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Abstract 

A  prism  chromotomographic  hyperspectral  imaging  sensor  is  being  developed  to 
aid  in  the  study  of  bomb  phenomenology.  Reliable  chromotomographic  reconstruction 
depends  on  accurate  knowledge  of  the  sensor  specific  point  spread  function  over  all 
wavelengths  of  interest.  The  purpose  of  this  research  is  to  generate  the  required  point 
spread  functions  using  wave  optics  techniques  and  a  phase  screen  model  of  system 
aberrations. 

The  phase  screens  are  generated  using  the  Richardson-Lucy  algorithm  for 
extracting  point  spread  functions  and  Gerchberg-Saxton  algorithm  for  phase  retrieval. 
These  phase  screens  are  verified  by  comparing  the  modeled  results  of  a  blackbody  source 
with  measurements  made  using  a  chromotomographic  sensor.  The  sensor  itself  is 
constructed  as  part  of  this  research.  Comparison  between  the  measured  and  simulated 
results  is  based  upon  the  noise  statistics  of  the  measured  image. 

Four  comparisons  between  measured  and  modeled  data,  each  made  at  a  different 
prism  rotation  angle,  provide  the  basis  for  the  conclusions  of  this  research.  Based  on 
these  results,  the  phase  screen  technique  appears  to  be  valid  so  long  as  constraints  are 
placed  on  the  field  of  view  and  spectral  region  over  which  the  screens  are  applied. 
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DESIGN  AND  MODEL  VERIFICATION  OF  AN  INFRARED 


CHROMOTOMOGRAPHIC  IMAGING  SYSTEM 


I.  Introduction 


Motivation 

Bomb  classification  using  remote  sensing  techniques  presents  a  significant 
research  challenge.  Recent  research  at  AFIT  has  linked  the  evolution  of  detonation 
spectra  in  time  to  the  initial  state  of  the  ordinance  with  some  success.  (Orson  and  others, 
2003: 107)  In  the  future,  these  classification  efforts  using  non- imaging  techniques  may 
be  supplemented  by  analysis  using  hyperspectral  techniques,  specifically 
chromotomographic  (CT)  imaging  technology. 

The  effectiveness  of  the  CT  imaging  system,  in  terms  of  the  spectral  and  spatial 
quality  of  the  data  product,  is  ultimately  tied  to  the  effectiveness  of  the  data 
deconvolution  algorithm.  The  original  algorithm  presented  in  the  literature  is  incapable 
of  reconstructing  spectral  information  from  images  dominated  by  low  spatial  frequency 
structure  (Mooney,  1997:2954).  To  improve  on  this  method,  a  new  deconvolution 
algorithm  is  being  developed  which  requires  an  explicit  understanding  of  the  point  spread 
function  (PSF)  specific  to  the  sensor. 
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Purpose 

The  purpose  of  this  research  is  to  assemble  a  CT  imager,  extract  the  phase 
information  required  to  build  appropriate  point  spread  functions,  and  demonstrate  the 
validity  of  these  phase  models  through  simulation  and  comparison  with  measured  CT 
imagery. 

Scope 

The  mid- wave  infrared  CT  sensor  is  assembled  and  configured  for  use  primarily 
in  the  laboratory  to  provide  a  validation  of  the  wave  optics  simulation  rather  than  a  fully 
deployable  sensor.  Similarly,  the  phase  retrieval  algorithm  used  to  generalize  the  PSF  is 
intended  to  demonstrate  the  concept  rather  than  provide  a  referendum  on  phase  retrieval 
algorithms.  The  simulation  approximates  the  propagation  of  light  in  tenns  of  Fourier 
transforms. 

Organization 

Progression  of  this  document  follows,  for  all  intents  and  purposes,  the 
chronological  development  of  the  research.  Chapter  II  contains  background  information 
on  hyperspectral  imaging  systems  and  specific  background  on  the  CT  system.  The  CT 
background  is  broken  down  into  principles  of  chromotomographic  imaging  and  a 
historical  prospective  of  CT  development. 

Chapter  III  explains  the  sensor  design  and  construction  process.  Prism  selection 
and  general  optical  system  criteria  are  explained  in  terms  of  geometric  optics.  This 
chapter  also  contains  an  explanation  of  the  component  selection  process  and  a  detailed 
description  of  sensor  assembly. 
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Chapter  IV  contains  the  theory  behind  the  propagation  model.  The  chapter  begins 
with  an  explanation  of  Fourier  propagation  and  image  formation  in  terms  of  the  point 
spread  function.  Later,  the  role  of  the  aberration  phase  screen  is  described.  The  chapter 
concludes  with  an  explanation  of  the  prism  transformation  model  and  composite  image 
formation.  Other  topics  of  interest  include  spatial  and  spectral  sampling. 

Chapter  V  deals  with  the  two  algorithms  required  to  generate  phase  screens  from 
measurements  of  intensity.  This  discussion  includes  the  laboratory  setup  required  to 
extract  a  point  spread  function  using  the  Richardson-Lucy  algorithm  as  well  as  a 
description  of  the  algorithm  itself.  The  second  algorithm,  the  Gerchberg-Saxton  phase 
retrieval  algorithm,  is  discussed  in  terms  of  its  functional  fonn  and  how  it  generalizes  the 
point  spread  function  extracted  using  the  Richardson-Lucy  algorithm. 

Chapter  VI  describes  the  noise  statistics  of  the  camera,  the  data  measurement 
process,  and  the  parameters  used  to  fit  the  data.  The  statistical  argument  is  based  upon 
the  assertions  of  other  researchers  as  well  as  experimental  data  specific  to  the  camera. 

The  measurement  process  is  presented  together  with  the  fit  parameters  because  the 
methodology  behind  one  helps  to  explain  the  other. 

Chapter  VII  contains  the  experimental  results  with  analysis  and  the  conclusions 
reached  during  this  research.  Four  comparisons  of  measured  and  simulated  broadband 
CT  images  are  used  in  this  analysis.  The  results  of  each  comparison  are  expressed  in 
terms  of  the  statistics  of  the  measured  image.  Further  discussion  is  devoted  to  the 
parameters  used  in  the  CT  model  to  approximate  the  measured  data  in  each  of  the  four 
cases. 
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II.  Background 


Hyperspectral  Imaging  Systems 

Traditional  non-imaging  spectrometers  provide  a  spectrum  formed  as  a 
conglomeration  of  emission  from  all  points  in  the  sensor’s  field  of  view.  A  Fourier 
transform  interferometer,  such  as  the  one  used  in  the  bomb  phenomenology  research 
referenced  previously,  is  an  example  of  this  type  of  sensor.  By  contrast,  hyperspectral 
imaging  (HSI)  systems  perform  the  same  function  over  the  field  of  view  of  each 
individual  pixel  in  the  sensor’s  focal  plane  array.  The  result  is  that  each  pixel  in  an  HSI 
formed  image  contains  its  own  independent  spectrum.  Hyperspectral  and  multispectral 
sensors  come  in  two  common  varieties,  those  in  the  filter  and  scanning  slit  category,  and 
a  third  emerging  variety,  the  chromotomographic  imager. 

The  product  provided  by  the  HSI  sensor  at  any  instant  in  time  is  commonly 
referred  to  as  the  data  cube.  Two  dimensions  of  this  cube  form  a  plane  that  corresponds 
to  the  spatial  dimensions  of  the  sensor’s  overall  field  of  view  and,  consequently,  each 
pixel  in  this  plane  has  a  field  of  view  equal  to  the  total  sensor  field  of  view  divided  by  the 
total  number  of  pixels  (assuming  100%  fill  factor).  The  value  assigned  to  each  pixel  in 
this  plane  represents  the  flux  incident  on  the  sensor  over  a  small  portion  of  the  spectrum. 
The  cube’s  third  dimension  identifies  spectral  position  and  has  a  length  equal  to  the 
spectral  bandwidth  of  the  sensor  divided  by  the  spectral  resolution  of  the  sensor.  Figure  1 
is  a  graphical  representation  of  the  data  cube  concept. 
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Figure  1:  The  I  lypcrspcctral  Data  Cube  (Modified  from  ITRES,2004) 


The  cube  concept  is  useful  for  describing  how  hyperspectral  data  can  be 
exploited.  When  the  cube  is  integrated  along  the  X  axis,  the  resulting  image  is  akin  to  a 
traditional  black  and  white  photograph.  Alternatively,  the  response  of  one  pixel  can  be 
plotted  against  X  and  the  result  is  a  spectrum  similar  to  what  would  be  obtained  by  a  non¬ 
imaging  spectrometer.  Furthermore,  individual  horizontal  slices  or  integrated  slices  can 
be  searched  separately  for  identifying  features  that  are  not  readily  apparent  in  traditional 
imagery.  The  potential  value  of  these  types  of  analysis  is  compounded  when  the  scene  in 
question  is  allowed  to  evolve  in  time  though  the  collection  of  a  series  of  data  cubes. 
Principles  of  the  Direct  Vision  Prism  CT 

A  polychromatic,  collimated  signal  enters  the  prism  assembly  and  is  dispersed 
about  a  central,  undeviated  wavelength.  The  dispersed  signal  is  collected  via  a  focusing 
lens  and  projected  upon  a  focal  plane  array.  This  process  is  updated  continuously  as  the 
prisms  rotate  but  the  sensor  is  otherwise  held  fixed  on  the  scene.  As  the  prisms  turns, 
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individual  frames  of  the  dispersed  image  are  collected  and  stored  for  analysis. 
Essentially,  each  recorded  frame  contains  a  superposition  of  many  monochromatic  slices 
of  the  image  offset  from  their  original  position  based  upon  the  dispersive  properties  and 
angle  of  rotation  of  the  direct  vision  prism.  Figure  2  is  a  graphical  representation  of  this 
process. 


rotatabt* 
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local  plana  array 
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•  peclral 
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Image  to 
computer 


by  turning  the  prism  we  get  a  sequence  ol  Images 


computer  is  used  to  separata  the  spectral  bands 


Figure  2.  CT  Image  Formation  Process  (Mooney,  1995:66) 


Once  the  image  sequence  is  collected,  the  convoluted  image  can  be  rebuilt  using 
an  image  reconstruction  algorithm.  Though  reconstruction  algorithms  vary  in  the  details, 
most  rely  on  the  notion  that  the  dispersed  image  can  be  expressed  mathematically  as  a 
“linear  superposition  of  each  spectral  image  convolved  with  a  unique  point  spread 
function.”  (Mooney,  1995:  67)  Equipped  then  with  the  dispersion  and  rotation  angles  of 
the  prisms  and  specific  knowledge  of  the  point  spread  function,  the  hyperspectral  data 
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cube  can  be  reconstructed  despite  the  unavoidable  presence  of  noise,  which  is  a 
significant  complicating  factor. 

Historical  Perspective  on  the  CT 

Though  a  variety  of  CT  imaging  systems  were  developed  somewhat  earlier,  the 
direct  vision  prism  MWIR  chromotomographic  camera,  which  is  the  basis  for  this 
research,  was  first  developed  by  Jonathan  Mooney  at  Rome  Laboratory,  Hanscom  Air 
Force  Base,  in  the  mid  1990’s.  (Mooney  and  others,  1997:2951)  Mooney  states  that  his 
development  trades  computational  overhead  (the  image  reconstruction  algorithm  requires 
significant  resources)  for  relatively  high  photon  throughput  and  temporal  resolution  when 
compared  to  modern  filter  or  slit  based  hyperspectral  imaging  systems.  Mooney’s 
reported  design  separated  the  3-5  um  portion  of  the  MWIR  band  into  25  sub-bands  over  a 
100  x  100  pixel  focal  plane  array. 

More  recently,  James  Murguia  of  Solid  State  Scientific  Corporation  (in 
conjunction  with  Mooney  and  others)  published  a  description  of  a  visible/near  IR  CT 
imager  capable  of  producing  64  spectral  bands  with  a  frame  rate  of  10  Hz.  (Murguia, 
2000:457)  Currently,  Sold  State  Scientific  commercially  offers  visible,  MWIR,  and 
LWIR  versions  of  the  CT  via  their  website.  (CTHIS,  2004) 

Two  AFIT  theses  on  CT  imaging  precede  this  work.  The  first  was  an  attempt  to 
simulate  basic  point  spread  functions  generated  by  a  notional  chromotomographic  imager 
to  aid  in  sensor  design  and  image  reconstruction  (these  are  goals  not  unlike  those 
presented  here,  but  on  a  much  more  modest  scale).  (Dearinger,  2004)  The  second  work 
was  a  characterization  of  two  candidate  reconstruction  algorithms  with  specific  emphasis 
on  the  recovery  of  absolute  radiometric  data.  (Gustke,  2004) 
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III.  CT  Design  Parameters  and  Selection 


Anatomy  of  the  Sensor 

Figure  3  contains  a  schematic  of  the  CT.  Dimensions  labeled  Dx  represent  the 
distances  between  objects  and  dimensions  Fx  represent  lens  focal  lengths. 


Figure  3.  The  CT  System  (Including  Dimensions) 

Moving  from  the  source  to  the  focal  plane,  a  brief  discussion  of  the  components 
of  the  CT  is  in  order  before  the  parameters  that  govern  them  are  described.  The  telescope 
has  three  primary  purposes:  to  collimate  the  input,  regulate  image  magnification,  and 
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ensure  that  the  signal  is  not  interrupted  by  the  edges  of  the  direct  vision  prism  system. 
After  collimation,  the  signal  enters  the  prism  assembly  where  it  is  dispersed  about  some 
center  wavelength  according  to  its  spectral  content.  During  operation,  the  prism 
assembly  rotates  and  so  does  the  dispersed  signal.  As  its  name  suggests,  the  focusing 
lens  recreates  the  dispersed  image  of  the  source  on  the  focal  plane. 

Prism  Development 

Constraints  on  the  Direct  Vision  Prism. 

Mr.  Kevin  Gross,  a  PhD  candidate  at  AFIT,  is  the  author  of  the  original  CT  prism 
design  in  a  collection  of  unpublished  work.  The  following  discussion  is  an  expansion  of 
his  work,  which  included  the  development  of  prism  design  constraints  and  a  survey  of 
appropriate  materials.  Based  on  his  conclusions,  the  fore  prism  (leftmost  in  figure  3)  is 
best  cut  from  Lithium  Floride  and  the  aft  prism  from  Barium  Floride.  His  original  design 
called  for  a  center  wavelength  at  3.6  lum,  the  center  of  the  InSb  bandpass,  but  this 
decision  was  made  before  it  was  known  that  the  camera  has  a  filter  that  narrows  the 
bandpass  down  to  3-5  urn.  The  proceeding  discussion  follows  the  development  of  the 
4um  center  wavelength  prism  set  designed  to  closely  match  the  attributes  of  the  original 
design.  Aside  from  the  enumeration  of  constraints,  what  follows  is  a  completely 
independent  development  from  Mr.  Gross’s  original. 

DVP Design  Constrains 

1)  assembly  must  fit  in  a  2  inch  diameter  circular  mount 

2)  undeviated  wavelength  should  be  4um 

3)  undeviated  wavelength  should  be  tolerant  of  cut  and  alignment  errors 

4)  dispersion  should  be  monotonic  over  the  3-5um  bandpass 

5)  useful  aperture  should  exceed  25%  of  prism  area 

6)  attenuation  and  reflection  losses  should  be  minimized 

7)  resolving  power  should  be  linear,  or  nearly  so 
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Center  Wavelength  and  Dispersion. 

Understanding  the  behavior  of  radiation  in  the  prism  system  is  a  prerequisite  for 
understanding  both  the  design  criteria  and,  later,  the  wave  optics  model  of  the  prism. 
Sellmeier  index  of  refraction  models  for  Barium  Floride  and  Lithium  Floride  used 
throughout  this  research  are  taken  from  (Troph,  1995: 1369).  Figure  4  is  a  diagram  of  the 
useful  quantities  in  the  following  derivation. 


Figure  4.  DVP  Ray  Trace  Diagram 


The  smaller  representation  of  the  system  identifies  the  relative  location  of  the  cutaway.  9i 
through  08  correspond  to  angles  used  throughout  the  derivation,  a  and  P  represent  prism 
cut  angles,  and  ni  through  n3  represent  the  media  specific  index  of  refraction. 

Radiation  entering  the  DVP  is  collimated.  The  radiation  path  through  the  system 
is  best  expressed  in  terms  of  successive  iterations  of  Snell’s  law. 


0,  =sin 


n, 


A 


-sin  a 


Vn2 


(3.1) 
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(3.2) 


02  =a-0j 


0,  =  sin  1 


n 


—  sin0. 


Vn3  J 

05=P-04 


07  =  sin  1 


— sin05^ 
vni  J 


08  =07  -P 


(3-3) 

(3.4) 

(3.5) 

(3.6) 


This  result  is  valid  with  or  without  the  presence  of  a  gap  between  the  prisms.  If  the  gap  is 
included 
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(3.8) 


hence  equation  3.17  is  identical  to  equation  3.12  and  the  rest  of  the  derivation  from  3.12 
on  continues  unabated. 

Sensitivity  to  Alignment  Error. 

The  preceding  development  assumes  that  the  two  prism  system  is  in  perfect 
alignment  though  realistically  this  may  not  always  be  the  case.  Net  alignment  error  can 
be  expressed  in  terms  of  the  effect  along  the  axis  of  dispersion  by  modifying  equations 
(3.4)  and  (3.6). 

05=Pcos(p-04  (3.9) 
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08  =  07  -PcOS(p 


(3.10) 


where  cp  is  the  angular  deviation  from  perfect  alignment  in  terms  of  the  aft  prism. 
Misalignment  also  gives  rise  to  dispersion  along  the  axis  perpendicular  to  the  axis  of 
intended  dispersion.  This  off-axis  dispersion,08p,  can  be  expressed  as 


08p  =  shT1 


—  sin  (P  sin  (p) 


(3.11) 


If  08p  is  ignored,  misalignment  can  be  interpreted  as  a  shift  in  center  wavelength. 
Useable  Aperture. 

The  left  of  figure  5  graphically  depicts  the  regions  1  through  5  referred  to  in  the 
following  development.  Figure  5  on  the  right  contains  the  end  on  view  of  the  prisms. 


1-2  Lj  L4 

II  II 


Figure  5.  Regions  of  the  Direct  Vision  Prism 

The  variable  r(x)  represents  the  total  prism  height  and  v(x)  is  the  maximum  ray 
entrance  height  possible  to  avoid  contact  with  the  prism  base.  Li  through  L3  represent 
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region  widths.  The  prisms  are  circular  therefore  the  values  r  and  v,  as  depicted  in  the 
cross  section  shown  in  figure  5  are  functions  of  x.  r(x)  is  defined  as 


r(x)  =  rmax  sin 


f  x  ' 

7T - 

V  max  J 


(3.12) 


Dimension  x  is  defined  as  distance  from  the  left  center  of  the  prism  cross-section. 

The  total  drop,  5,  from  entrance  to  exit  is  expressed  as  the  sum  of  the  drops  in 
each  region. 

5  =  (r(x)-  v(x))-  v(x)tanatan02  -L1  tan02  -L2  tan0g  -L3  tan04  (3.13) 


The  value  of  v  is  detennined  by  setting  5  equal  to  zero  and  solving.  Inspection  of 
equation  3.17  reveals  that  5  is  independent  of  P  but  not  of  n3.  In  this  context,  percent 
useful  aperture,  UA,  is  defined  as 


UA  = 


|  v(x) • dx 
| r(x) • dx 


(3.14) 


Note  that  the  integral  over  r(x)  is  simply  the  surface  area  of  the  prism  face.  The  angles 
specified  in  equation  3.19  are  detennined  at  a  given  wavelength.  This  wavelength  should 
be  selected  to  coincide  with  maximum  dispersion  if  UA  is  to  represent  the  region  where 
all  radiation  across  the  bandpass  is  permitted  to  pass  through  the  entire  prism  system 
without  reflection. 

Attenuation  and  Reflection  Losses. 

The  radiation  entering  the  prism  is  unpolarized.  Transmittance  across  each 
interface  can  be  expressed  in  terms  of  the  Fresnel  equations  (Hecht,  2002: 1 15-120) 


13 


Tr±(0i,0t,ni,nt)  = 


^nt  cos0t  ' 

2sin0t  cos  0j 

ln;  cos0,J 

v  sin(0j  +0t)  y 

(3.15) 


Tr=(0i,0t,ni,nt)  = 


cos0t  ' 

2sin0tcos0i 

Ui  cos0,J 

vsin(0i+0t)cos(01-0t)y 

(3.16) 


In  each  case,  0i  and  0t  represent  incident  and  transmitted  angles  and  n,  -  nx  represent  the 
respective  index  of  refraction.  The  subscripts  on  Tr  represent  the  orientation  of 
polarization  (parallel  and  perpendicular).  Total  transmission  through  the  prism  system  is 
given  as 


Tr  (total)  = 


(3.17) 


The  index  4  represents  the  four  interfaces  encountered  during  propagation. 

Resolving  Power. 

The  sensor’s  ability  to  discriminate  spectral  features  is  an  important  attribute  to 
understand.  One  possible  measure  of  this  ability  is  resolving  power  (Pedrotti,  1987: 122) 


R 


X 

(AT.) 


(3.18) 


where  X  represents  wavelength  and  the  denominator  represents  the  minimum  discemable 
distance  (in  terms  of  wavelength)  between  spectral  features.  Pedrotti’s  complete 
development  is  particular  to  a  single  prism  specified  by  its  apex  angle  and  is  therefore  not 
appropriate  for  the  direct  vision  prism  though  what  follows  is  a  development  in  the  same 
spirit.  The  Rayleigh  distance,  the  distance  between  the  peak  of  the  diffraction-limited 
point  spread  function  and  the  first  null  serves  as  a  gage  for  detennination  of  R.  The 
Rayleigh  range  (Pedrotti,  1987:335)  is  defined  as 
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(3.19) 


DR(A) 


1.22- A -F3 
h 


Note  that  each  variable  is  defined  as  before.  Using  DR,  (AA)m;n  can  be  expressed  in  the 
following  terms 


Dr  (A,  +  ais-)  +  Dr  (A.  -  sa-)  =  1  ~22 '  2A 

R  2  R  2  h 


DR  (A  +  AAmin )  +  DR  (A  -  AAmin )  < 


F3  tan  08  (A  +  -  F3  tan  08  (A  - 


2.44A 

h 


08(A  +  ^l)-0s(A-^ 


(3.20) 

(3.21) 

(3.22) 


In  this  case,  08  is  expressed  as  a  function  of  A  and  the  small  angle  approximation 
for  tangents  is  applied  in  equation  3.2 1 .  (AA)min  is  best  detennined  numerically  by 
representing  08(A)  as  a  third  order  polynomial. 

There  are  two  limits  to  the  resolving  power  model.  First,  this  definition  ceases  to 
have  meaning  if  the  point  spread  function  falls  entirely  in  one  pixel  on  the  focal  plane 
array.  Second,  the  results  (as  stated)  are  accurate  only  in  the  diffraction-limited  case  and 
hence  this  result  is  an  absolute  upper  limit  on  resolving  power. 

Attributes  of  the  Selected  Prisms 

Center  Wavelength  and  Dispersion. 

Figure  6  contains  the  direct  vision  prism  dispersion  angles  with  respect  to 
wavelength  for  both  the  3.61  um  and  4.00  um  center  wavelength  designs.  The  width  of 
each  curve  represents  the  error  associated  with  the  manufacturer’s  prism  cut  tolerances 
(±30  arcsec).  Dispersion  over  the  3  to  5  um  band  is  not  linear  in  either  case  but  is 
monotonic.  The  total  dispersion  about  the  bandpass  is  less  than  1  degree.  The  4  um 


15 


curve  is  well  described  as  a  linear  shift  in  wavelength  from  the  3.61  um  curve,  a  property 
which  manifests  itself  when  calculating  resolving  power. 


Figure  6.  Direct  Vision  Prism  Dispersion  Angles 
Curves  representing  a  ±5°  alignment  error  stand  out  to  the  right  of  the  perfect  alignment 
curves  for  each  set  of  prisms.  Though  not  shown  on  the  plot,  this  error  also  produces  an 
approximately  constant  shift  of  0.03°  along  an  axis  perpendicular  to  the  primary 
dispersion. 

Useful  Aperture. 

Both  prism  systems  exhibit  a  useful  aperture  of  92%,  a  result  guaranteed  by 
choosing  the  same  material  and  cut  angle  for  the  leading  prism.  Figure  7  represents  the 
concept  of  useful  aperture  graphically.  In  the  plot,  5  um  radiation,  the  wavelength  of 
largest  dispersion,  is  reflected  off  of  the  base  of  the  prism  assembly  when  it  falls  in  the 
region  between  the  r(x)  and  v(x)  curves.  This  region  of  reflection  forms  a  crescent  along 
the  edge  of  the  leading  prism  at  the  base. 
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x,  the  distance  along  the  prism  center  (mm) 

Figure  7.  Graphical  Representation  of  the  Useful  Aperture 

Attenuation  and  Reflection  Losses. 

Total  transmission  through  prism  varies  between  0.884  at  3  um  to  0.900  at  5  um 
for  both  prism  designs.  Like  the  dispersion  function,  the  transmission  function  is 
monotonic  and  otherwise  unremarkable. 

Resolving  Power. 

The  plot  in  figure  8  shows  that  resolving  power  increases  linearly  with  increasing 
wavelength.  By  inspection  of  equation  (3.18),  this  result  is  an  indication  of  constant 
spectral  resolution  over  the  prescribed  bandpass.  Though  only  one  curve  is  plotted 
below,  the  resolving  power  results  for  both  prism  sets  are  virtually  indistinguishable  (as 
predicted  in  the  previous  discussion  on  dispersion). 


—  r(x) 

-  V(x) 
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Figure  8.  Resolving  Power  as  a  Function  of  Wavelength 


Lens  Development 

Constraints  on  Lens  Selection. 

In  general,  lens  selection  constraints  are  restricted  only  what  is  required  for 
functionality.  The  fourth  and  fifth  items  in  the  following  list  of  constraints  are 
constraints  determined  both  by  lens  selection  and  by  prism  selection  (refer  ahead  to 
equation  3.8). 


Lens  Selection  Constraints 

1)  lenses  are  selected  from  an  array  of  2  inch  CaF  lenses  specified  at  f/l  through  f/10 

2)  collimated  light  must  pass  through  the  DVP  without  interacting  with  the  mount  walls 

3)  field  of  view  must  be  large  enough  to  image  a  1  inch  diameter  source 

4)  theoretical  propagation  theory  predicted  field  of  view  should  be  maximized 

5)  image  magnification  must  not  result  in  signal  falling  off  the  focal  plane 

6)  required  bandwidth  must  fall  entirely  on  the  focal  plane 
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Matrix  Representation  of  the  Optical  System. 

Ignoring  the  prisms  for  the  moment,  the  transformation  matrix  of  the  system  can 
be  expressed  as  such  using  the  following  treatment.  (Pedrotti,  1987:  70)  Translations 
take  the  fonn 


T(DX) 


(3.23) 


and  the  transformation  at  a  lens  interface  is  defined  as 


1  0 


L(FX)  = 


1 


(3.24) 


In  terms  of  these  two  components,  the  transformation  of  the  system,  M,  can  be  defined  as 
such 


M  =  T(F3 )  •  L(F3 )  •  T(D3 )  •  L(F2 )  •  T(D2 )  •  L(F1 )  •  T(D1 )  (3.25) 

The  system  defined  by  M  will  form  an  image  on  the  focal  plane  when  the  following 
condition  is  satisfied 


M12 


D,(D2-F1-F2)  +  F1(F2-D2)  0 

f,f2 


(3.26) 


where  the  subscripts  on  M  represent  the  row  and  column  position  of  the  term  in  question. 
Equation  3.4  can  be  solved  for  D2 


D 


2 


F>iF1+D1F2-F1F2 

Di-F, 


(3.27) 
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Notice  that,  in  the  limit  where  Di  goes  to  infinity,  Do  =  F1+F2,  which  is  the  expected 
result.  Finally,  we  can  establish  the  transverse  magnification  of  the  system  by 
substituting  this  result  for  Do  into  Mu 


Mu  = 


_p  Ft  +  F2  -  P2 


(3.28) 


Both  the  expressions  for  magnification  and  D2  are  completely  independent  of  the  distance 
D3  and  hence  this  dimension  is  discarded  as  a  design  constraint. 

Field  of  View. 

The  field  of  view,  0fov  is  an  expression  of  the  angular  width  of  the  scene  captured 
on  the  FPA.  Using  similar  triangles,  the  FOV  can  be  expressed  in  terms  of  the  linear 
width  of  the  focal  plane,  W,  and  the  effective  focal  length  of  the  system 


(3.29) 


If  the  entire  optical  system  is  modeled  as  a  single  thin  lens,  the  effective  focal  length,  feff, 
is  given  by  manipulating  the  well  known  thin  lens  equation 


1 


- h 

D, 


1 


(3.30) 


Note  that  the  effective  focal  length  defined  here  is  not  the  effective  focal  length 
associated  with  principal  plane  analysis,  it  is  simply  the  calculated  focal  length  of  the 
system  if  it  were  modeled  as  a  single  thin  lens  located  at  the  real  position  of  the  leading 
optic  in  the  CT  system.  For  now,  the  importance  of  feff  can  be  summed  up  via  the 
following  expression  for  the  Fresnel  propagation  region  of  validity  (Goodman,  1996:69) 


(feff)3»^[(X-tf  +(y-fi)2]  (3-31) 
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where  (x,y)  and  (2,r|)  represent  the  position  in  the  image  and  object  planes  respectively. 
As  later  discussion  will  show,  (3.3 1)  may  be  overly  restrictive  but  the  advantage  of 
maximizing  feff  is  clear. 

Determination  of  the  Focal  Lengths  Ft  and  F2. 

Ideally,  the  signal,  once  it  has  passed  through  the  telescope,  will  have  no 
interaction  with  the  walls  of  the  prism  mount.  The  radiation  from  the  source,  under  all 
conceivable  operating  conditions,  can  be  considered  nearly  collimated  at  the  aperture  of 
the  sensor.  Under  these  conditions,  the  telescope  is  approximately  afocal  and  the 
relationship  between  the  diameter  of  the  objective  lens,  H,  and  the  diameter  of  its  image 
(the  exit  pupil),  h,  can  be  approximated  as  follows  (Hecht,  2002:221) 


H  =  F1 

h  F2 


(3.32) 


Assuming  the  conditions  on  the  field  of  view  are  met,  the  problem  of  selecting  Fi 
and  F2  is  therefore  reduced  to  ensuring  three  things:  h  is  less  than  the  diameter  of  the 
useful  prism  aperture,  the  resulting  system  magnification  is  sustainable,  and  the  effective 
focal  length  is  as  large  as  possible.  Note  that  the  last  two  requirements  represent 
competing  interests.  Stray  radiance  from  around  the  edges  of  the  objective  lens  is 
prevented  from  entering  the  system  via  a  stop  placed  between  the  lenses  in  the  telescope. 
In  this  configuration,  this  stop  serves  as  the  field  stop  of  the  system. 
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Factors  Influencing  the  Selection  of  the  Focusing  Lens. 

The  focusing  lens  defines  both  the  magnification  of  the  system  from  equation  3.6 
and  the  linear  dispersion  of  the  image  which,  depending  on  the  dimensions  of  the  focal 
plane,  can  limit  bandwidth.  Given  a  prism  exit  angle,  0g,  the  linear  dispersion,  L,  at  the 
focal  plane  is  simply 

L  =  F3  tan 08  (3.33) 


Implicit  in  this  expression  is  the  counterintuitive  requirement  that  neither  D3  nor  the 
distance  between  the  prism  and  the  focusing  lens  affects  the  linear  dispersion  at  the  focal 
plane. 

Under  certain  circumstances,  pixel  size  will  also  effect  the  selection  of  F3. 

Though  it  is  clear  from  (3.22)  that  resolving  power  is  not  a  function  of  F3  in  an  analog 
sense,  (3.19)  shows  that  the  width  of  the  diffraction  limited  PSF  is  tied  to  F3.  In  the  limit 
where  PSF’s  become  small  compared  to  pixel  width,  the  definition  of  resolving  power 
losses  validity.  This  concern  should  be  recognized  in  general  but  can  otherwise  be 
largely  ignored  because,  as  the  results  of  this  research  will  show,  even  the  smallest  PSF’s 
have  a  width  measurable  over  several  pixels. 

Transmission  Losses  in  the  Lenses. 

Losses  in  the  lenses  are  assumed  to  occur  due  to  reflections  at  normal  incidence  at 
each  of  the  six  air/glass  or  glass/air  interfaces.  Transmittance  at  nonnal  incidence  is 
given  by  (Hecht,  2002:121) 


Tr  =  4"'n' 


n  /  \2 

(nt+n.) 


(3.34) 
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where  nt  and  n,  represent  the  index  of  refraction  of  the  incident  and  transmitted  media. 
Since  it  is  clear  than  transmittance  at  normal  incidence  is  equivalent  at  both  air/glass  and 
glass/air  interfaces,  total  transmittance  through  three  lenses  is  given  by  the  sixth  power  of 
Trn.  These  losses  are  insubstantial,  unavoidable  and  are  not  necessary  as  part  of  the 
design  process. 

Attributes  of  the  Selected  Lenses 

Telescope  Lens  Selection. 

As  stated  in  the  constraints,  lens  selection  must  be  placed  in  the  context  of  the 
selected  prisms.  The  direct  vision  prism  system  has  a  diameter  of  25  mm  and  a  useful 
aperture  of  approximately  91%.  By  equation  (3.32),  the  ideal  ratio  of  Fi  to  F2  would  be 
slightly  less  than  2. 1  to  fully  utilize  the  useful  aperture  though  prism  interior  reflections 
will  result  if  this  ratio  is  reduced  further.  A  ratio  of  2.5  is  attainable  using  the  f/3  and  f/2 
lens  combination  as  is  ratio  of  3  using  an  f/3  and  f/1  combination.  Of  the  two  choices, 
effective  focal  length  receives  a  boost  by  choosing  the  f/3  -  f/1  combination. 
Consequently,  the  f/3  (15.24  cm  focal  length)  and  f/1  (5.08  cm  focal  length)  combination 
is  selected  to  be  the  best  compromise. 

Focusing  Lens  Selection. 

There  are  four  primary  factors  that  influence  the  selection  of  F3:  bandwidth, 
magnification,  effective  focal  length  and  field  of  view.  The  results  of  the  angular 
dispersion  plot  in  figure  6  combined  with  equation  (3.33)  shows  that  the  total  linear 
dispersion  between  3  to  5  um  on  the  focal  plane  will  cover  228  pixels  (each  pixel  is  24 
um  wide)  using  the  f/10  lens.  The  entire  focal  plane  has  dimensions  of  640  by  5 12 
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pixels  therefore,  even  in  the  most  extreme  case,  selection  of  focal  length  F3  will  not 
restrict  the  desired  bandwidth. 

Effective  focal  length  and  total  magnification  both  increase  with  an  increase  in 
F3.  The  case  for  increasing  feff  has  already  been  made  but  field  of  view  must  be  kept  in 
check  to  prevent  portions  of  the  incident  signal  critical  to  PSF  extraction  from  leaking  off 
the  focal  plane.  The  source  described  in  lens  selection  criteria  only  requires  a  field  of 
view  of  0.33°  if  it  were  imaged  perfectly  but  aberrations  expand  this  requirement 
significantly.  Based  on  the  comparison  in  table  2  and  with  this  restriction  in  mind,  the  f/5 
lens  (25.4  cm  focal  length)  is  best  suited  for  this  research. 

Summary  of  Selection  Results 

Tables  3  and  4  summarize  the  selection  of  optical  components  and  their  attributes. 
Aside  from  center  wavelength,  the  attributes  of  both  prism  designs  are  virtually 
indistinguishable.  Though  interesting,  this  result  is  not  necessarily  surprising  because  the 
center  wavelength  change  was  implemented  by  only  a  small  change  in  the  Barium 
Floride  prism  cut  angle.  For  comparison,  several  results  for  the  f/10  focusing  lens  are 
also  included. 


Table  1.  Specifications  of  the  Direct  Vision  Prism 


Center 

A  (um) 

angle  a 

angle  p 

Total 

Dispersion 

UA 

Transmission 

(%) 

Resolving 

Power 

3.61 

18°33'06" 

14°34'30" 

-0.16°  to  0.46° 

0.92 

0.88  to  0.90 

40  to  73  (linear) 

4.00 

18°33'06" 

14°20'34" 

-0.26°  to  0.35° 

0.92 

0.88  to  0.90 

40  to  73(linear) 
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The  total  dispersion  defined  in  table  1  is  in  terms  of  the  3  to  5  um  bandpass.  The  two 
prism  systems  are,  in  terms  of  most  figures  of  merit,  identical. 


Table  2.  Lens  System  Specifications 


Fi  (cm) 

F2  (cm) 

F3  (cm) 

Dt  (cm) 

D2  (cm) 

Magnification 

F  effective 
(cm) 

FOV 

(deg) 

15.24 

5.08 

25.40 

434.34 

20.87 

0.18 

66.82 

0.81 

15.24 

5.08 

50.80 

434.34 

20.87 

0.36 

115.82 

0.55 

The  first  row  of  table  2  contains  the  specifications  of  lenses  used  during  this  research,  the 
second  row  exists  for  the  sake  of  comparison. 

Assembly  and  Alignment  of  the  Direct  Vision  Prism  System 
DVP  Assembly. 

Material  considerations  complicate  what  would  otherwise  be  the  simple  task  of 
coupling  and  mounting  the  direct  vision  prisms.  According  to  the  prisms’  manufacturer, 
optically  bonding  the  Barium  Floride  and  Lithium  Floride  prisms  is  possible,  but  no 
adhesive  is  available  with  favorable  transmission  properties  in  the  mid-wave  1R. 
Furthermore,  holding  the  two  prisms  together  with  a  heat  activated  shrink  wrap  is  also 
impossible  because  Barium  Floride  is  extremely  sensitive  to  thermal  shock.  As  a  result, 
the  two  prisms  are  separately  mounted  mechanically  and  then  attached  and  aligned. 
Figure  9  shows  the  prisms  in  their  mounts,  coupled  together,  and  attached  to  the  rotation 
stage. 
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Figure  9.  The  Direct  Vision  Prism  Assembly 


The  fore  prism  mount  (Lithium  Floride)  consists  of  a  cylindrical  section  of 
machined  aluminum  with  a  concentric  inner  ring  cut  to  hold  the  prism.  The  prism  is 
wrapped  in  an  adhesive  backed  neoprene  sheath  and  held  against  a  lip  facing  the  aft 
prism  section.  The  sheath  is  necessary  to  prevent  chipping  at  the  glass-aluminum 
interface  (a  lesson  learned).  A  set  screw  is  used  to  secure  the  wrapped  prism  in  place. 
This  mount  is  held  to  the  aft  prism  mount  via  a  bolt  in  a  hole  that  sweeps  out  a  10°  arc  for 
alignment  purposes. 

The  aft  prism  mount  also  contains  a  circular  hole  and  set  screw  system  to  contain 
the  prism  and  hold  it  against  a  lip  facing  the  fore  prism  mount.  This  mount  bolts  directly 
to  the  rotation  stage.  To  aid  the  explanation,  figure  10  contains  a  diagram  of  the  two 
mounts. 
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Figure  10.  Schematic  of  the  Prism  Mounts 

The  two  lips  create  a  net  air  gap  between  the  two  prisms  of  0. 16  cm.  The  presence  of  the 
air  gap  prevents  the  two  prisms  from  rubbing  during  alignment  but  also  reduces  the 
useful  aperture  (the  calculations  in  table  1  include  this  gap). 

D  VP  Alignment 

Inspection  of  the  geometry  in  figure  3  indicates  that,  when  the  prisms  are 
perfectly  aligned,  both  reflections  and  transmissions  through  the  prisms  will  occur  in  the 
plane  whose  nonnal  is  perpendicular  to  the  optical  axis.  Analytically,  this  statement  is 
reinforced  by  the  development  for  cp  in  the  previous  discussion  of  alignment  errors.  A 
NeHe  pointing  laser,  shown  in  figure  9,  is  used  to  exploit  this  relationship  for  alignment 
purposes.  Proper  alignment  is  achieved  when  two  points  of  reflection,  captured  on  a 
paper  screen  facing  the  fore  prism  assembly,  and  two  points  of  transmission,  captured  on 
a  screen  facing  the  aft  assembly,  all  fall  in  the  expected  plane.  Alignment  deviations 
show  up  as  a  rotation  out  of  this  plane  away  from  the  optic  axis  nonnal.  Figure  9  shows 
the  bright  primary  transmission  point  and  a  dim  secondary  reflection  dot  lined  up  on  the 
aft  prism  facing  screen. 
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IV.  The  Fourier  Transform  Propagation  Model 


Each  simulated  broadband  image  is  formed  as  a  summation  of  many 
monochromatic  images.  In  turn,  these  monochromatic  images  are  formed  as  a 
convolution  of  the  geometrically  predicted  image  of  the  source  with  a  point  spread 
function  specific  to  the  sensor.  The  following  development  explains  how  these  point 
spread  functions  are  derived  and  applied  for  the  ultimate  purpose  of  image  fonnation. 
Propagation  Theory  Fundamentals 
Fresnel  Propagation. 

The  following  expression  relates  the  field  in  the  image  plane  to  the  field  in  the 
object  plane  using  the  Fresnel  propagation  integral.  (Goodman,  1995:67) 


U(x,  y)  =  —  e 

j^Z 


j27iz 

E(x2+y2)  “  * 


A.z'' 
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uftn)< 


S^) 


*A,z 


d^dp 


(4.1) 


where  U(x,y)  and  U(^,p)  represent  the  field  in  the  image  and  object  planes  and,  for 
purposes  of  simplicity,  each  field  has  units  of  the  square  root  of  intensity  (v  Watts  /  cm)  . 
The  variable  z  in  (4. 1 )  represents  the  distance  between  the  planes  described  by  (x,y)  and 
(^,q).  By  applying  the  following  change  of  variables: 


fx  =  — 

A,z 

f  =  *- 
y  Xz 


(4.2) 


the  Fresnel  integral  can  be  transfonned  into  a  scaled  Fourier  transform 


U(x,  y)  = 


j27TZ 


jXz 


(4.3) 
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where  the  script  F  represents  the  Fourier  transform  operator.  In  both  cases,  the 
corresponding  intensity  pattern  is  given  by 

I(x,y)  =  U(x,y)U*(x,y)  (4.4) 

Multiplication  by  the  conjugate  annihilates  the  effects  of  any  constant  phase  terms. 
Consequently,  it  is  often  convenient  to  drop  these  terms  early  on  during  computation  to 
simplify  bookkeeping. 

Two  Essential  Properties  of  the  Fresnel  Propagation. 

A  limit  on  the  validity  of  Fresnel  propagation  was  imposed  earlier  in  equation 
(3.3 1)  in  terms  of  the  effective  focal  length.  In  general,  this  equation  can  be  recast  in 
tenns  of  z  but,  as  explained  in  the  reference  (Goodman,  1996:69),  this  limit  does  not 
apply  in  all  cases.  As  a  result,  quantifying  the  region  in  the  image  plane  where  the 
Fresnel  approximation  applies  in  difficult  to  determine  in  the  absence  of  system  specific 
measurements.  Any  number  of  separate  Fresnel  regions  could  be  defined  for  the  focal 
plane  of  the  CT  system  but  this  is  a  task  best  left  to  a  future  thesis.  Assuming  the  images 
made  by  the  CT  fall  into  a  single  one  of  these  regions,  the  following  statements  define  a 
system  a  properties  shared  by  every  included  intensity  pattern.  Except  where  specifically 
noted,  reference  for  these  properties  can  be  found  in  chapter  2  of  the  Goodman  text. 

Invariance. 

An  optical  system  is  considered  to  be  temporally  and  spatially  invariant  if  the 
only  result  of  a  change  in  an  object’s  position  is  a  change  in  its  image.  Note  that  this 
definition  applies  specifically  to  the  optical  system  but  makes  no  reference  to  a  source, 
which  may  or  may  not  be  invariant.  This  is  a  property  shared  by  both  Fresnel  and 
geometric  propagation.  Put  another  way,  if  I(x,y)  is  the  image  formed  from  an  object 
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described  by  I(^,q)  and  a  is  a  proportionality  constant  that  links  one  coordinate  system  to 
the  other  then  for  a  space  invariant  system 

Iferi)i»8ft-g0,ri-ri0)  >I(x,y)®8(x-^0,x-an0)  (4.5) 

where  20  and  r|0  represent  the  position  of  the  center  of  the  shifted  object  intensity  pattern 
and  5  symbolizes  the  Dirac  delta  function.  Though  the  reasons  why  are  not  explicitly  laid 
out  in  equation  (4.5),  this  property  will  be  critical  to  the  development  of  the  prism 
transformation  phase  screen. 

Linearity. 

A  generic  linear  operator,  S,  exhibits  the  following  properties  in  all  cases 

S{u-I1(x,y)  +  v-I2(x,y)}  =u-S{l1(x,y)}  + v-S{l2(x,y)}  (4.6) 

where  u  and  v  are  scalar  modifiers  of  intensity  patterns  fi  and  I2.  The  operator  S  could 
represent  the  Fresnel  propagation  itself  (which  leads  directly  to  the  following  property)  or 
S  could  represent  a  pattern  summing  operator  which  allows  for  the  superposition  of 
intensity  patterns  to  form  a  single  image.  Consequently,  this  property  can  be  used  to 
describe  a  polychromic  intensity  pattern  as  a  composite  of  many  individual 
monochromatic  patterns. 

Image  Formation  as  a  Convolution 

The  image  formed  by  a  linear  invariant  optical  system  can  be  expressed  as  a 
convolution  of  the  system  specific  impulse  response  (point  spread  function)  and  the 
geometrically  predicted  image  of  the  source.  (Goodman,  1995:  21) 

I(x,  y)  =  s(x,  y)  0  g(x,  y)  (4.7) 
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where  s(x,y)  is  the  normalized  intensity  pattern  of  the  point  spread  function  and  g(x,y)  is 
the  geometric  intensity  pattern  of  the  source. 

Calculation  of  g(x,y)  is  straightforward  if  the  origins  of  both  image  plane  (x,y) 
and  the  object  plane  (^,q)  are  defined  to  lie  on  the  optical  axis. 

g(x,y)  =  g(-|— ,-7— )  (4.8) 

Mu  Mn 

where  Mu  is  the  magnification  defined  in  (3.28).  This  expression  has  nothing  to  do  with 
Fresnel  propagation;  it  is  simply  a  mapping  from  one  plane  to  the  other. 

Calculating  the  Point  Spread  Function 

Physically,  the  point  spread  function  is  formed  as  the  image  of  a  unit  intensity 
point  object.  Analytically,  this  field  corresponding  to  the  point  object  at  the  aperture, 
Ap(2,r|),  that  has  the  following  property 

00  00 

J  J  Ap(§,  r|)  •  Ap(^,  r|)*d^dr|  =  1  (4.9) 

-00  -00 

Equation  (4.9)  states  that  the  aperture  must  have  a  finite  area  but  places  no  restrictions  on 
the  dimensions  of  that  area.  Furthermore,  this  definition  places  no  restrictions  on  the 
phase  at  the  aperture  (which  includes  atmospheric  aberrations  induced  during 
propagation).  Regardless  of  the  size  of  the  aperture,  the  total  intensity  collected  is 
nonnalized  to  1 . 

A  lens  is  required  to  form  the  image.  The  function  form  of  a  thin  lens  is  taken  to 
be  (Goodman,  1995:99) 

-kt-(^2+r12) 

t1(4,fi)  =  e  Xf  (4.10) 
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where  f  is  the  effective  focal  length  of  the  lens.  This  phase  model  is  adequate  to  form  a 
theoretical  image  but  real  images  include  lens  aberrations  which  requires  the  addition  of 
another  phase  tenn  (Roggemann,  1996:27) 

tab(^)  =  e*  (4.11) 


where  y(^,q)  is  a  measurement  of  the  phase  aberrations  at  every  point  in  the  pupil, 
presented  in  terms  of  wavelengths.  This  aberration  function  can  be  build  as  linear 
combination  of  individual  aberrations  using  Zemike  polynomials  (Roggemann,  1996:95) 
or  they  can  be  sampled  from  the  optical  system  itself  (a  process  which  is  demonstrated  in 
this  research).  tab  can  also  be  used  to  represent  near  field  atmospheric  aberrations  that 
occur  during  propagation  through  space.  For  brevity  sake,  all  lens  transformations  will 
henceforth  be  expressed  as  a  combination  of  the  thin  lens  and  aberration  transfonnations. 


t(M  =  e 


kf 


(if+r> 


(4.12) 


All  the  tools  required  to  form  the  point  spread  function  are  now  in  place. 

Referring  back  to  (4.3)  the  field  produced  in  the  focal  plane  due  to  the  point  object  is 
given  to  be 

j27TZ 

UPSF(x,y)  =  ^— eJ^fx  +fy  ^(25^,|Ap(4,'n)t(4,'n)e^  +11  ^  j  (4.13) 

Completely  expanded,  two  of  the  exponential  terms  cancel  because,  for  focused  image 
formation  z  must  equal  f,  and  the  equation  reduces  to 


UPSF(x,  y) 


e  x  jjiXz(fx2+fv2)  r  \  ,  ,  ^r(5.n)l 

—  e  U  Ca  j AP(^,r|)ex 


(4.14) 
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and  hence  the  point  spread  function  intensity  pattern  is  given  by 

s(x,  y)  =  UPSF(x,  y)UPSF*(x,  y)  (4.15) 

A  Single  Lens  Phase  Model 

At  first  glance,  it  appears  that  at  least  four  Fresnel  propagations  are  necessary  to 
model  the  CT  optical  system  since  there  are  three  separate  lenses  between  the  source  and 
the  focal  plane.  A  multi-propagation  approach  will  yield  the  correct  result  but  the  system 
is  greatly  simplified  by  the  constraint  that  radiation  entering  the  direct  vision  prism  is 
collimated  (or  nearly  so).  This  constraint  is  manifest  in  equation  (3.25)  where  the 
distance  between  lens  F3  is  set  to  be  the  length  F3. 

If  the  edge  effects  of  the  optics  in  the  interval  are  approximately  nil,  the  source 
can  be  considered  to  be  an  aberrated  plane  wave  incident  on  the  focusing  lens.  These 
aberrations  are  picked  up  in  the  atmosphere,  in  the  telescope  portion  of  the  optical 
system,  and  in  the  prisms,  but,  as  long  as  they  are  represented  adequately  (sampled 
properly),  they  need  not  be  considered  separately.  Of  course,  the  efficacy  of  this 
approximation  is  dependent  upon  the  acceptability  of  its  effect  on  the  results. 

As  stated  previously,  this  approximation  is  not  absolutely  necessary  but  it  is 
preferable  for  two  reasons.  The  first  (and  lesser)  reason  is  computational  simplicity;  the 
number  of  Fourier  Transforms  required  for  a  given  propagation  is  reduced  from  one  for 
each  optical  interface  to  a  total  of  one.  Second,  while  it  is  possible  to  measure  the  effects 
of  aberrations  separately  in  each  optical  element,  unless  perfect  alignment  is  achieved, 
the  net  effect  of  these  aberrations  is  not  equal  to  the  sum  of  its  parts.  Treating  the  all 
aberrations  as  an  event  in  a  single  plane  and  making  corresponding  measurements 
accounts  for  both  optical  imperfections  and  alignment  errors  while  approximating  out  the 
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edge  effects  of  the  intervening  optics.  Making  measurements  of  these  aberrations  is  the 
subject  of  the  proceeding  chapter. 

Definition  of  the  Aperture  Function. 

The  field  associated  with  collimated  light  is  described  analytically  in  terms  of  a 
plane  wave.  Furthermore,  the  collimated  beam  at  the  exit  of  the  telescope  has  a  diameter 
h  as  defined  in  equation  (3.32).  In  terms  of  plane  waves  alone,  the  normalized  aperture 
function  is  defined  as 


APpla„e(X3>y3) 


j27tD3 


Vx7~+y7  -  h 


tt(^) 

0  ^x32  +  y32  >  h 


(4.16) 


where  (X3,y3)  are  the  coordinates  in  the  plane  of  the  lens  F3.  Note  that  the  constant  phase 
tenn  in  this  equation  does  depend  on  D3  but,  as  stated  previously,  constant  phase  is  of  no 
consequence  in  terms  of  intensity.  The  value  of  h  corresponds  to  the  radius  of  the 
aperture. 

Aberrations  introduced  into  the  signal  due  to  atmosphere  and  imperfect  telescope 
optics  are  lumped  together  as  a  virtual  thin  lens  applied  to  the  system  in  the  plane  of 
propagation.  These  aberrations  combined  with  the  previous  result  form  the  field  at  the 
entrance  of  the  focusing  lens 


Ap(x3,y3)  =  Applane 


(Scx3,Scy3)e 


V3) 


(4.17) 


where  yprop  is  the  atmosphere  and  telescope  aberration  phase  screen  and  Sc  is  a 
wavelength  scaling  factor  that  will  be  described  in  the  next  section. 
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Physical  Manifestation  of  yprop. 

A  brief  digression  is  useful  at  this  point  to  provide  context  for  equation  (4.17).  If 
YproP(x3,y3)  were  set  to  zero  at  all  points  and  the  subsequent  aperture  function  substituted 
into  the  transformation  in  (4. 14),  the  resulting  field  at  the  focal  plane  (see  figure  11) 
would  be  the  familiar  Airy  pattern,  the  diffraction  limited  point  spread  function. 


Figure  11.  The  Diffraction  Limited  PSF  (Hecht,  2002:469) 

In  cases  where  yproP(x3,y3)  is  not  zero  over  all  space,  the  convolution  property  of 
the  Fourier  transform  is  applied  (Goodman,  1996:9). 

C/lAp^e’’”  1  =  C/{APP,„„}  ®  c/pS'’-'|  (4. 18) 

Hence  the  field  in  the  focal  plane  can  be  interpreted  as  the  convolution  of  the  Airy  disk 
with  the  transform  of  the  aberrations  accumulated  along  the  optical  path. 

The  phase  screen  Yprop  is  measured  at  a  single  wavelength  at  a  given  aperture  size. 
The  magnitude  of  the  aberration  is  achromatic  therefore  a  scaling  factor  is  introduced  into 
the  exponent  so  that  the  phase  measurements  can  be  applied  to  all  wavelengths  of 
interest. 
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(4.19) 


where  X  is  the  wavelength  of  propagation  and  Am  is  the  wavelength  where  the  phase 
screen  was  initially  measured.  The  similarity  theorem  of  the  Fourier  transform 

C/jg(Scx,Sty)}=TGf|L>|Lj  (4.20) 

demonstrates  that  the  effect  of  this  scaling,  besides  nonnalizing  the  phase  screen,  is  to 
scale  the  result  of  the  transform  in  image  space  by  the  inverse  of  the  scale  factor.  The 
physical  implication  of  this  result  is  that  the  PSF  will  be  wider  for  longer  wavelengths 
and  more  narrow  for  shorter  wavelengths.  This  contraction  and  expansion  is  verified 
independently  by  the  definition  of  the  Rayleigh  Criterion  in  equation  (3.19). 

The  Direct  Vision  Prism 

Up  to  this  point,  no  consideration  has  been  given  to  the  effects  of  the  direct  vision 
prism  on  the  model.  Recall  from  equation  (4.5)  the  requirement  for  spatial  invariance  in 
the  image  plane.  Under  this  regime,  the  effect  of  the  prism  must  be  to  shift  the  image 
(and  by  logical  extension,  the  PSF)  but  leave  it  otherwise  completely  intact.  The  effect 
on  the  PSF  due  to  the  prism  in  the  image  plane  is 

U’psf  =UpSF®8(x-u,y-v)  (4.21) 


where  u  and  v  represent  linear  dispersion  in  the  image  plane.  The  variables  u  and  v  share 
the  following  relationship  with  equation  (3.33) 


u  =  F3  tan  08p 
v  =  F3  tan  08 


(4.22) 
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The  effects  of  dispersion  are  brought  back  into  the  diffraction  plane  (2,r|)  by  the  inverse 
Fourier  transform  of  (4.21). 


C/'/{U’pSF}  =  Applanee 


fScY, 


prop  A^(5tan08p+ritaneg) 


(4.23) 


In  tenns  of  phase,  the  prism  is  manifests  itself  as  an  additional  aberration  term  in  the 
plane  of  diffraction. 

The  Discrete  Direct  Vision  Prism  Transformation. 

The  preceding  derivation  expresses  all  transformations  in  terms  of  continuous 
functions.  With  the  exception  of  the  prism  transformation,  the  continuous  case  migrates 
into  the  discrete  case  with  relative  ease,  assuming  that  specific  sampling  requirements  are 
met  (see  the  following  section  on  sampling).  The  discrete  prism  transformation  is 
analogous  to  its  continuous  counterpart  but  sufficiently  different  the  two  warrant  some 
specific  attention. 

The  circular  shift  property  of  the  DFT  (Strum,  1989:399)  is  expressed  in  one 
dimension  as 


j2n 


D]g(k)eN  l  =  G(n©m) 


(4.24) 


where  D  is  the  DFT  operator,  g(n)  is  the  transformation  of  G(k),  N  is  the  number  of 
samples  in  the  array,  and  m  is  the  specified  shift.  For  this  application,  m  corresponds  to 
the  direction  and  magnitude  of  the  linear  dispersion  in  a  particular  direction 

m  =  —  (off  -  axis) 

dx  (4.25) 

m  =  —  (on -axis) 
dx 
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where  u  and  v  come  from  the  development  in  (4.22)  and  dx  is  the  sampling  distance  in 
the  image  plane.  Consequently,  equation  (4.23)  can  be  recast  in  discrete  space  as 

jn  j2n( k,F3  tanQgp  |k2F3tan9sN| 

D'jUVSF(n„n2)}=Ap[l„1,(k1,k,)e^'v'eNl  *  a’  1  (4.26) 

where  the  indices  nx  and  kx  represent  the  sampled  versions  of  the  corresponding 
functions.  At  this  point,  the  parallel  between  (4.23)  and  (4.26)  is  obvious. 

Formation  of  the  Composite  Image  and  Radiometry 

Equation  (4.26)  rounds  out  all  of  the  required  steps  necessary  to  generate  a  point 
spread  function  for  any  wavelength  of  interest.  Additionally,  the  discrete  version  of 
equation  (4.7)  provides  the  mechanism  necessary  to  apply  the  PSF  to  the  problem  of  real 
image  fonnation.  What  remains  then  is  to  combine  these  results  to  form  a  complete 
polychromatic  image  by  appropriately  scaling  the  intensity  in  terms  of  wavelength  and 
summing  the  resulting  intensity  patterns. 

Scaling  with  Intensity. 

In  terms  of  the  geometrically  predicted  image,  the  total  power  per  unit  wavelength 
incident  on  the  focal  plane,  P(k),  is  given  by 

p(^)=ZZgfeii)A^Aii  (4-27) 

where  A^Aq  is  the  area  of  a  pixel  in  the  focal  plane  and  g  is  defined  as  in  equation  (4.8). 
This  equation  provides  a  simple  mechanism  for  ensuring  that  the  correct  amount  of 
power  is  deposited  on  each  pixel  in  tenns  of  each  wavelength  of  interest.  Throughout 
this  research  the  source  is  a  circular  cavity  blackbody  of  temperature  T.  The  problem 
remains  then  to  define  power  per  unit  wavelength  incident  on  the  focal  plane: 
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P(X)  = 


Lbb(k,T)AsAo 


(4.28) 


f 

V 
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TfTcTwTaTCaF 


where  the  tenn  in  parentheses  is  a  measure  of  the  total  power  from  a  source  (of  area  As) 
collected  by  the  leading  optic  of  the  system  (area  A0)  over  a  distance  of  Di.  The  tenn  Tx 
represent  the  various  mediums  and  filters  though  which  the  radiation  must  pass.  In  order 
of  appearance  from  left  to  right  these  absorbers  include: 


Tf  a  source  filter  (set  to  1  for  broadband  sources) 

Tc  transmission  through  the  camera  cold  filter 

Tw  transmission  through  the  camera  window 

Ta  transmission  through  the  atmosphere  (Modtran) 

TCaF  transmission  through  the  three  Calcium  Floride  Lenses 


Each  of  these  transmission  terms  deserves  some  discussion  but  first  all  qualms 
about  lumping  these  tenns  together  must  be  laid  to  rest.  Recall  from  the  discussion  of 
linearity  of  the  Fourier  transform  that  constant  tenns,  in  this  case  the  magnitude  of  the 
electric  field  or  one  of  its  modifiers,  may  be  pulled  out  of  each  transformation  and 
considered  separately.  As  a  result,  the  magnitude  electric  field  at  the  image  plane  (and 
hence  the  total  power)  need  only  to  be  calculated  once  as  the  product  of  the  propagated 
field  and  its  modifiers,  which,  in  this  case  are  the  various  absorbers. 

All  sources  suffer  the  effects  of  the  last  five  absorbers  mentioned  above.  Data  on 
Tc  and  Tw  is  taken  from  the  camera  manufacturers  documentation  (SBFP,  unk:  13-54)  and 
transmission  through  the  Calcium  Floride  lenses  comes  from  the  treatment  in  equation 
(3.34)  along  with  the  subsequent  discussion. 

Approximation  of  the  transmission  through  the  atmosphere  was  determined  using 
Modtran4’s  standard  US  atmospheric  model  at  constant  pressure  with  a  path  length  of  5 
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meters.  The  actual  path  length  from  source  to  FPA  is  slightly  less  than  this  distance  but 
Modtran4  does  not  support  path  length  increments  of  less  than  1  meter.  For  more 
information  on  Modtran4  see  (Berk,  2003). 

Rounding  out  equation  (4.28),  Lbb(k,T)  represents  blackbody  photon  radiance  of 
the  source  (Dereniak,  1996:66) 


where 


for  the 


Lbb(^T)  =  - 


2c 


f  he  "\ 


(4.29) 


k4 


em-  _  ] 


h  is  Planck’s  constant,  c  is  the  speed  of  light,  and  k  is  Boltzman’s  constant. 

A  typical  example  of  P(k)  for  a  blackbody  source  at  400°C  is  given  in  figure  12 
3  to  5  um  bandpass. 
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Figure  12.  Blackbody  Spectral  Flux  at  the  Focal  Plane 


The  upper  (smooth)  curve  represents  the  total  photon  flux  in  the  absence  of  absorbers  and  is 
provided  as  a  reference.  The  lower  curve  is  the  flux  used  to  approximate  P(A.)  in  the  simulation. 
Notice  the  prominent  atmospheric  feature  found  near  4.3  um.  The  edges  of  the  bandpass  a 
defined  primarily  by  the  cold  stop  filter. 
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Summing  Intensity  Patterns. 

The  second  issue  is  also  easily  resolved.  By  definition  this  optical  system  is 
linear  in  intensity.  As  a  result,  a  polychromatic  image  can  be  expressed  as  a  sum  of  many 
(approximately)  monochromatic  images  assuming  that  the  image  is  sampled 
appropriately  in  wavelength.  The  Nyquist  sampling  theorem  (Strum,  1988:54) 
guarantees  that  the  analog  polychromatic  signal  can  be  recovered  so  long  as  the  sampling 
frequency,  fs,  meets  the  following  criterion: 

fs  >  2fmax  (4.30) 

where  fmax  is  the  zero  amplitude  cutoff  frequency  of  the  power  spectrum.  Note  that  the 
power  spectrum  is  the  Fourier  transform  of  P(k).  In  practice,  a  zero  amplitude  frequency 
may  not  exist  but  can  be  approximated  by  selecting  the  cutoff  where  the  amplitude  is 
several  orders  of  magnitude  lower  that  the  peak.  Figure  13  contains  a  plot  of  the  Fourier 
transform  of  P(k)  for  both  a  simple  blackbody  at  400°C  (the  smooth  curve)  and  the 
blackbody  modified  by  the  aforementioned  absorbers.  As  a  reference,  10  urn-1  can  be 
considered  an  adequate  cutoff  for  sampling  purposes. 
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Figure  13.  Fourier  Transform  of  the  Blackbody  Spectrum 


Given  this  cutoff,  fs  must  be  at  least  20  unT1  and  hence  sampling  should  occur  at  0.05  um 
intervals. 

Prism  Rotation 

The  rotation  of  the  direct  vision  prism  can  be  handled  in  one  of  two  ways 
depending  on  the  geometry  of  the  source.  For  circularly  symmetric  sources,  such  as  the 
blackbody  source  used  in  this  research,  prism  rotation  can  be  applied  directly  to  the 
composite  intensity  pattern,  I(x,y),  using  the  following  transformation  matrix 


x’~ 

e°s0rot 

sin0rot 

X 

y'_ 

-sin0rot 

cos  0rot  _ 

_y_ 

where  9rot  is  the  prism  rotation  angle  measured  counterclockwise  from  the  line  x  =  0. 

Most  sources  in  nature  will  not  be  circularly  symmetric.  In  these  cases,  the 
transformation  in  (4.31)  should  be  applied  directly  to  the  prism  phase  screen  before  it  is 
multiplied  with  the  aberration  phase  screen.  This  approach  is  also  suitable  for  circular 
sources  but  much  less  efficient. 
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V.  Phase  Screen  Calculation 


The  previous  section  discusses  the  application  of  ypr0p  but  the  problem  of  deriving 
Ypr0p  remains.  First,  a  concise  explanation  of  what  the  phase  screen  does  and  does  not 
contain  is  appropriate.  Recall  that  yprop  describes  the  diffraction  plane  field  phase  of  the 
on-axis  space  and  time  invariant  point  spread  function.  Consequently,  yprop  is 
independent  of  the  source.  Additionally,  no  mention  has  been  given  to  the  field  effects  of 
the  prism  assembly;  the  screen  is  intended  to  capture  only  the  effects  of  aberrations 
caused  by  lenses  and  their  alignment.  Given  these  restrictions,  the  first  task  is  to  extract 
the  point  spread  function  from  measurements  made  in  the  laboratory. 

True  monochromatic  point  sources,  which  are  required  for  a  direct  measurement 
of  the  point  spread  function,  are  difficult  to  fabricate  in  the  laboratory.  To  circumvent 
this  problem,  the  Richardson-Lucy  algorithm  is  applied  to  a  measurement  of  a  (nearly) 
monochromatic  source  with  known  spatial  dimensions  to  statistically  estimate  a  PSF  at 
the  center  wavelength.  To  generalize  the  PSF  to  all  wavelengths,  a  second  algorithm  is 
then  applied  to  extract  the  phase. 

The  Richardson-Lucy  Algorithm 

The  Richardson-Lucy  algorithm  (Lanteri,  1972:55)  addresses  the  problem  of 
recovering  an  original  image  from  its  degraded  measurement.  In  Lanteri’s  parlance,  the 
relationship  between  the  two  is  given  as 

f(x,y)  =  g(x,y)®h(x,y)  (5.1) 

where  f  is  the  degraded  image,  g  is  the  original  image,  and  h  is  the  point  spread  function. 
Consequently,  the  problem  of  recovering  g  amounts  to  a  deconvolution.  Given  f(x,y)  and 
h(x,y),  the  Richardson-Lucy  algorithm  suggests  an  iterative  approach  to  this  problem: 
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(5.2) 


gk+i  (x,  y)  =  gk  (*,  y)  • h*  (-x,  -y)  ® 


f(x,y) 

h(x,y)®gk(x,y) 


where  k  is  the  index  of  iteration  and  go(x,y),  the  initial  guess,  is  constant,  positive,  and 
nonnalized. 

Without  appropriate  modification,  the  Richardson-Lucy  algorithm  does  not 
provide  a  mechanism  for  calculating  h(x,y),  which  is  the  stated  goal.  As  an  intennediate 
step,  the  commutative  property  of  convolution  must  be  invoked 

f(x,y)  =  h(x,y)®g(x,y)  (5.3) 

and  the  Richardson-Lucy  algorithm  is  recast  as 

ht+,(x.y)=M*.y)-g'(-*.-y)«>  ,  f)*’y),  ,  (5.4) 

g(x,y)®hk(x,y) 

which  is  an  iterative  solution  for  estimating  the  point  spread  function  given  the  original 
image  and  its  corresponding  measurement.  As  an  aside,  note  that  equation  (5.3)  is 
identical  to  equation  (4.7). 

Source  Setup. 

The  RL  algorithm  for  calculating  h(x,y)  requires  as  input  both  the  perfect  (ray 
traced)  image  of  a  source,  g(x,y),  and  a  measurement  of  the  monochromatic  source  image 
though  the  optics,  f(x,y).  A  monochromatic  source  is  approximated  in  the  lab  using  a 
blackbody  coupled  with  a  narrowband  thin-film  filter.  Figure  14  contains  the 
measurement  of  this  particular  filter’s  response  made  using  a  Bomem  MR-254  FTIR 
spectrometer.  Radiation  passing  through  the  filter  will  fall  predominantly  between  4.2 
and  4.4  um. 
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Figure  14.  Filtered  Blackbody  Irradiance,  Normalized  to  the  Unfiltered  Response 

Sampling  Issues. 

Like  the  wavelength  sampling  issues  discussed  in  chapter  4,  similar  issues  must 
also  be  addressed  here  in  two  dimensions.  The  ultimate  goal  is  to  adequately  sample  the 
point  spread  function  which,  by  extension,  requires  that  both  g(x,y)  and  f(x,y)  are 
sampled  correctly.  The  pixel  size  in  the  focal  plane  determines  the  sampling  in  f(x,y)  and 
is  therefore  also  a  convenient  sample  size  to  choose  for  g(x,y).  To  test  if  this  sample  size 
is  adequate,  the  Fourier  transfonn  of  both  g(x,y)  and  f(x,y)  are  examined  to  ensure  that  a 
reasonable  cutoff  frequency  is  reached.  Figure  15  contains  a  measurement  of  f(x,y)  (left) 
and  its  Fourier  Transform  (right). 
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Figure  15.  A  Captured  Image  f(x,y)  and  its  Fourier  Transform 


The  image  is  of  a  filtered,  circular  blackbody  source  of  diameter  1.02  cm  formed  as  the 
average  of  15  frames  taken  at  72  Hz.  The  background  (due  to  internal  reflections  and 
emission  along  the  optical  path)  is  suppressed  and  replaced  with  a  small  random 
background.  This  step  is  necessary  to  prevent  to  edge  of  the  camera  cold  stop  from  being 
interpreted  as  part  of  the  source  image.  Figure  16  contains  the  same  treatment  applied  to 
g(x,y). 
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Figure  16.  The  Generated  Image  g(x,y)  and  its  Fourier  Transform 
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The  leftmost  image  in  figure  16  contains  the  geometrically  predicted  image  of  the 
source  in  figure  15  and  the  rightmost  image  contains  its  Fourier  transform. 

The  two  image  functions  and  their  Fourier  transforms  provide  sufficient  evidence 
of  proper  sampling.  Both  transformed  functions  are  nonzero  over  only  a  finite  region  in 
Fourier  space.  Functions  that  exhibit  this  behavior  (or  closely  approximate  it)  are 
considered  to  be  bandlimited.  Functions  of  this  class  enjoy  a  special  property:  a  properly 
sampled  bandlimited  function  can  be  used  to  completely  reconstruct  the  original  analog 
function  (Goodman,  1996:  23).  In  general,  Fourier  transfonns  are  periodic.  Bandlimited 
functions  are  properly  sampled  when  this  periodic  structure  does  not  overlap  (a 
phenomenon  known  as  aliasing).  It  is  clear  from  figures  15  and  16  that  no  such  overlap 
occurs.  Consequently,  both  functions  are  adequately  sampled. 

Iterating  the  RL  Algorithm. 

The  Richardson-Lucy  algorithm  does  not  come  packaged  with  a  definitive  method 
for  detennining  the  number  of  required  iterations.  Qualitatively,  this  cutoff  should  occur 
once  subsequent  iterations  no  longer  improve  significantly  on  the  previous 
approximation.  To  quantify  this  concept,  the  sum  of  the  squared  error  between 
neighboring  approximations,  Esq,  is  introduced 

Esq  =  ZZ[hk(x’y)-hk+i(x5y)]2  (5-5) 

x  y 

Unless  Esq  goes  to  zero,  the  cutoff  iteration  is  left  to  some  interpretation.  Figure  17 
contains  a  semilog  plot  of  Esq  versus  iteration  number.  This  particular  data  comes  from 
the  application  of  the  Richardson-Lucy  algorithm  on  the  g(x,y)  and  f(x,y)  combination 
provided  in  the  previous  discussion. 
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Figure  17.  Sum  Squared  Error  of  the  Richardson-Lucy  Algorithm 


The  magnitude  of  the  squared  error,  large  or  small,  is  of  little  practical  use  (recall  the 
h(x,y)  is  normalized  and  is  therefore  unitless).  Instead,  the  cutoff  iteration  is  selected 
based  on  the  (somewhat  subjective)  determination  of  where  slope  of  the  Esq  curve  is 
effectively  zero.  At  500  iterations,  the  Esq  curve  presented  above  appears  to  meet  this 
criterion. 

The  Estimated  Point  Spread  Function. 

Figure  18  contains  the  image  of  the  PSF,  h5oo(x,y),  corresponding  to  the  work  laid 
out  in  the  previous  discussion. 
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Figure  18.  The  Richarchson-Lucy  Estimated  PSF 


Recall  that  this  is  an  estimation  of  the  PSF  at  around  4.3  um.  The  next  step  is  to 
create  a  generalized  phase  screen  out  of  this  information  that  can  be  used  to  describe 
point  spread  functions  at  all  wavelengths. 

The  Gerchberg-Saxton  Phase  Retrieval  Algorithm 

The  Gerchberg-Saxton  Phase  Retrieval  Algorithm  (Gerchberg,  197 1 :237) 
estimates  diffraction  plane  phase  given  measurements  of  the  intensity  in  the  diffraction 
and  imaging  planes.  For  this  specific  application,  the  point  spread  function  generated  by 
the  Richardson-Lucy  algorithm  provides  the  intensity  pattern  in  the  imaging  plane  and 
the  intensity  in  the  diffraction  plane  is  estimated  to  be  constant  with  respect  to  position. 
As  a  result,  the  sampling  scheme  used  to  implement  the  Richardson-Lucy  algorithm  is 
maintained. 

According  to  the  author,  the  GS  algorithm,  presented  schematically  in  figure  19, 
takes  advantage  of  the  fact  that  a  change  in  amplitude  in  one  plane  (diffraction  or  image) 
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is  the  result  of  a  change  in  both  phase  and  amplitude  in  the  other.  Gerchberg  goes  on  to 
show  that  the  squared  error  between  the  image  reconstructed  using  the  phase  from  his 
algorithm,  h|((x,y),  and  the  original  image  (in  this  case,  the  Richardson-Lucy  recovered 
PSF),  hri(x,y),  must  decrease  or  remain  constant  with  each  successive  iteration. 

Z  Z  [Mx>  y)  -  hk+i(x>  y)f  <=  Z  Z  [hri  (x>  y)  -  hk  (x>  y)]2  (5-6) 

x  y  x  y 

Note  that  the  subscript  k  has  been  recycled  to  again  refer  to  iteration  index. 


START 


Figure  19.  Schematic  Representation  of  the  Gerchberg- Saxton  Algorithm  (Gerchberg,  1972:239) 


Iterating  the  Gerchberg-Saxton  Algorithm. 

Equation  (5.6)  leaves  open  the  possibility  that  the  algorithm  will  stagnate  before 
perfect  reconstruction  of  the  point  spread  function  is  achieved  (i.e.  when  the  squared 
error  is  driven  to  zero).  Figure  20,  a  plot  of  this  squared  error  over  500  iterations,  shows 
that  the  Gerchberg-Saxon  reconstruction  of  the  PSF  hri(x,y)  is,  in  fact,  subject  to  this 
predicted  stagnation. 
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Figure  20.  Stagnation  of  the  Gerchberg- Saxton  Phase  Retrieval  Algorithm 


x  10 


Comparing  Recovered  Point  Spread  Functions 

Ultimately,  the  impact  of  this  imperfect  reconstruction  is  difficult  to  assess  out  of 
the  context  of  the  finished  product  but,  in  the  interval,  a  comparison  between  the  GS 
estimated  PSF  and  the  RL  recovered  PSF  will  suffice. 

Figure  21  contains  the  Richardson-Lucy  recovered  point  spread  function  from 
figure  18,  enlarged  and  cropped  so  that  a  direct  comparison  can  be  made  between  it  and 
the  Gerchberg-Saxton  estimated  point  spread  function. 
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Figure  21.  Enlarged  Richardson-Lucy  Point  Spread  Function 


For  comparison,  figure  22  contains  the  Gerchberg- Saxton  estimation  of  the  point  spread 
function  in  figure  2 1 . 
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Figure  22.  PSF  Generated  from  Gerchberg-Saxton  Recovered  Phase 


Visual  inspection  seems  to  indicate  a  strong  resemblance  between  the  two  figures 
but  the  two  are  clearly  not  identical.  Notably,  more  power  appears  to  be  concentrated  at 
the  center  of  the  image  in  figure  21  and,  by  comparison,  the  central  feature  in  figure  22  is 
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more  amorphous.  Predictably,  this  center  broadening  will  manifest  itself  in  the  CT  model 
as  a  more  dispersed  source  though,  as  stated  previously,  the  total  impact  of  this  problem 
is  best  left  for  the  final  analysis. 

Reiteration  the  entire  process  places  the  calculation  of  the  phase  screen  in  its 
proper  context.  The  recovered  phase,  yprop,  is  scaled  to  the  appropriate  wavelength  and 
combined  with  the  phase  generated  by  the  prism  model  to  produce  a  model  of  the  point 
spread  function  on  the  focal  plane  array.  The  mechanism  for  this  process  was  established 
in  chapter  IV.  The  geometric  image  of  the  source  is  convolved  with  the  point  spread 
function  to  produce  a  monochromatic  image  of  the  source.  Many  such  images  are  then 
combined  to  produce  a  complete  simulation  of  the  CT  response  to  the  source  input.  At 
this  point  the  model  is  complete  and  attention  can  turn  to  the  process  of  comparing 
simulated  images  with  imagery  made  in  the  laboratory. 


53 


VI.  Measurements  and  Mechanisms  for  Comparing  Results 


The  model  defined  in  chapters  IV  and  V,  though  in  general  quite  flexible,  places 
certain  requirements  on  how  laboratory  measurements  must  be  prepared  and  processed. 
These  requirements  fall  into  five  categories:  nonuniformity  correction,  background 
suppression,  image  alignment  errors,  prism  alignment  errors,  and  rotation  errors.  With 
the  exception  of  responsivity,  each  of  the  other  requirements  are  dealt  with  as  individual 
parameters.  Before  exploring  these  parameters,  a  sufficient  explanation  of  the  statistics 
of  an  image  captured  from  the  CT  is  in  order. 

Techniques  such  as  correlation  and  sum  squared  error  provide  a  means  for 
determining  which  simulation  parameters  best  fit  a  particular  result.  Both  of  these 
mechanisms  will  be  applied  in  this  analysis  but  neither  can  be  described  as  a  definitive 
test  for  measuring  the  overall  effectiveness  of  the  simulation.  An  excellent  example  of 
this  dilemma  can  be  made  by  reaching  ahead  into  chapter  VII  and  examining  some 
correlation  results.  Two  simulations  are  made  and  correlated  with  a  measurement 
resulting  in  correlation  coefficients  of  0.9133  and  0.9895.  Both  correlations  turn  out  to 
be  close  to  1  (1  being  complete  correlation)  but  one  of  the  simulated  images  has  been 
rotated  by  90°  from  the  image  it  was  intended  to  simulate.  The  two  results  appear 
favorable  because  the  intensity  pattern  changes  in  the  center  of  the  image  but  most  of  the 
pixels  considered  are  not  affected  by  the  rotation  and  hence  remain  highly  correlated. 
This  comparison  clearly  shows  that  the  90°  rotated  simulation  is  not  the  better  of  the  two 
choices  though  the  result  provides  no  insight  into  overall  agreement  with  the  collected 
data. 
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Image  Statistics 

Model  and  data  agreement  is  best  couched  in  terms  of  the  noise  statistics  of  the 
image.  The  values  recorded  at  each  pixel  are  a  quantized  representation  of  the  number  of 
incident  photons  being  absorbed  in  the  pixel’s  volume  over  some  interval.  The 
probability  of  the  camera’s  electronics  reporting  any  particular  value  is  governed  by 
either  photon,  quantization,  or  electronic  noise  depending  on  the  intensity  of  the  source 
and  the  quality  of  the  instrument.  Once  this  distribution  is  identified,  the  quality  of  the 
simulated  data  can  be  expressed,  pixel  by  pixel,  in  terms  of  the  probability  of  making  a 
measurement  that  exactly  matches  the  simulation.  The  outward  form  of  this  expression  is 
represented  by  the  number  of  standard  deviations  a  simulated  pixel  is  away  from  its 
measured  counterpart. 

The  Poisson  Distribution. 

According  to  E.  L.  Dereniak,  the  optical  system  described  above  should  follow 
Poisson  statistics: 

...photons  follow  Poisson  statistics  for  all  practical  detector  applications.  That  is, 
for  visible  and  near  infrared  applications,  low  temperature  blackbodies,  and  short 
wavelength  [radiation]  (hv  »  kT),  the  photon  noise  follows  Poisson  statistics 
(Dereniak,  1996:  156) 

As  a  result,  the  Poisson  distribution  is  used  as  a  theoretical  model  for  the  CT  camera 
image  statistics.  For  the  intended  purpose,  the  most  important  characteristic  of  this 
distribution  relates  the  mean,  p,  to  the  standard  deviation,  a. 

(Kreyszig,  1999:1081) 

p  =  a2  (5.7) 
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To  determine  the  applicability  of  this  distribution,  the  CT  optical  system  was 
configured  to  take  700  frames  (the  maximum  number  allowed  by  the  camera  buffer)  of  a 
temporally  constant  blackbody  source.  The  mean  and  standard  deviation  of  these  images 
is  measured  and  compared  to  a  theoretical  Poisson  distribution.  Figure  23  contains  a 
cross  section  of  the  focal  plane  array  though  the  center  of  the  source  image  along  with  the 
theoretical  and  measured  error  bars. 


4 
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Figure  23.  Image  Cross  Section  with  Theoretical  and  Measured  Error  Bars 

The  thick  solid  line  in  figure  23  is  the  mean  value  of  each  pixel,  enumerated  along 
the  x  axis  and  the  error  bars  are  identified  by  the  dashed  lines.  Short  dashes  identify  the 
theoretical  Poisson  error  bars  and  long  dashes  identify  the  actual  error  bars.  The  error  bar 
overlap  is  not  perfect,  but  given  the  sample  size  (700  frames)  the  difference  is  justifiable. 

Sufficient  evidence  of  a  Poisson  distributed  system  exists  between  the 
combination  of  the  collected  noise  data  presented  in  figure  23  and  Dereniak’s  assertions 
for  regions  where  signal  is  present.  Though  not  included  in  the  plot,  regions  of  the  focal 
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plane  where  little  or  no  signal  was  detected  did  not  follow  Poisson  statistics.  This  result 
is  not  surprising  because  other  fonns  of  noise  become  more  prevalent  in  the  absence  of 
signal  (i.e.  dark  current).  Failure  of  the  statistical  model  in  these  regions  is  acceptable 
because  they  have  little  impact  on  the  model  verification  process. 

Processing  Images 

The  noise  statistics  of  an  image  are  presented  first  so  that  the  effects  of  the  image 
processing  procedure  can  be  better  understood.  As  stated  previously,  the  model 
presented  in  this  research  requires  that  measurements  be  made  and  processed  in  a  certain 
way.  The  remainder  of  this  section  will  deal  with  this  processing  method  and,  where 
appropriate,  its  effects  on  the  statistical  fit  of  the  model. 

Nonuniformity  Correction. 

Pixel  response  is  a  combination  of  pixel  responsivity  and  the  connected 
amplification  electronics.  Nonuniformity  correction  (NUC)  is  applied  to  all  imaging 
systems  with  more  than  one  detector  element  to  compensate  for  variations  in  individual 
pixel  bias  and  gain.  This  process,  also  referred  to  as  flat  field  calibration,  is  carried  out 
by  exposing  the  entire  focal  plane  array  to  sources  of  uniform  intensity  and  adjusting  the 
gain  and  bias  of  each  pixel  so  that  their  response  is  equivalent  at  these  of  these  points. 
Figure  24  is  a  graphical  depiction  of  this  process. 
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Response  of  Pixel  2 


Figure  24.  Two  Pixel  Nonuniformity  Correction 


As  shown  in  the  graphic,  the  NUC  is  carried  out  by  matching  the  count  output  of 
two  pixels  at  two  known  levels  of  irradiance.  Nonuniformity  correction  is  applied  in  the 
CT  system  by  replacing  the  focusing  lens  with  a  unifonn  blackbody  source  and  making 
measurements  at  two  known  temperatures.  The  details  of  NUC  processing  are  carried  out 
by  software  provide  by  the  camera’s  manufacturer.  The  NUC  is  sensitive  to  thermal 
cycling  therefore  a  new  correction  is  made  each  time  the  camera  is  cooled  from  room  to 
liquid  nitrogen  temperature. 

Background  Suppression. 

The  model  makes  no  allowance  for  photons  generated  from  anywhere  besides  the 
intended  target  whereas,  in  the  laboratory,  it  is  impossible  to  isolate  the  intended  target  in 
a  single  measurement.  This  background  may  include  thennal  photons  generated  inside  of 
the  CT  system,  photons  scattered  onto  the  FPA  from  outside  the  system,  and  photons 
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from  the  scene  surrounding  the  target.  Eliminating  or  suppressing  these  background 
sources  is  dealt  with  though  a  combination  of  physical  and  software  means. 

The  effects  of  photons  generated  inside  of  the  CT  system  can  be  subtracted  away 
from  the  target  image  assuming  this  background  remains  constant  and  can  be  isolated 
from  the  scene.  The  internal  background  is  made  constant  by  shrouding  the  system  with 
black  velvet  cloth.  As  an  additional  precaution,  the  distance  between  the  focusing  lens 
and  the  camera  face  is  enveloped  by  a  cardboard  tube.  Furthennore,  the  region  in  front  of 
the  leading  aperture  is  enclosed  with  thick  packaging  material  to  eliminate  reflections 
from  the  laboratory  floor.  These  devices  emit  weakly  over  the  bandpass  but  the  effective 
internal  background  is  limited  to  only  to  this  emission.  Figure  25  is  a  picture  of  the 
shrouded  CT  system. 


Figure  25.  The  Shrouded  CT  System 

Having  isolated  the  sensor  from  its  surroundings,  the  next  step  is  to  isolate 
the  background  from  the  blackbody  source.  This  is  accomplished  taking  an  image  of  the 
scene  (with  the  source  target  covered)  seconds  before  the  target  image  is  made.  The 
background  image  contains  the  internally  and  externally  generated  background  but  not 
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the  source  or  its  immediate  surroundings  (the  blackbody  and  the  face  of  the  apparatus  that 
contains  it).  To  emphasize  the  importance  of  this  step,  figure  26  contains  a  slice  of  an 
image  made  before  background  subtraction  (left)  and  the  same  slice  of  an  image  after 
background  subtraction  (right). 
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Figure  26.  Images  Before  and  After  Background  Subtraction 


As  seen  in  figure  26,  a  small  bias  remains  after  background  subtraction.  This 
remaining  background  is  likely  the  result  of  thermal  emission  from  the  region 
immediately  surrounding  the  blackbody  cavity  that,  during  collection  of  the  background 
image,  is  obscured  along  with  the  cavity.  Suppression  of  this  effect  is  treated  as  an 
adjustable  parameter  whose  selection  is  based  upon  the  image  correlation  coefficient. 


C  = 


i=i  j=i 


XXM(i,j)2XXN(i,j)2 

f  i=l  j=l  i=l  j=l 


(5.8) 


where  M  and  N  represent  the  measured  and  simulated  images.  The  value  of  the 
suppression  parameter  is  selected  based  upon  its  ability  to  maximize  the  C.  The 
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correlation  coefficient  is  more  appropriate  for  selection  of  this  parameter  than  squared 
error  because  squared  error  emphasizes  intensity  differences  where  magnitudes  are  large 
while  the  correlation  coefficient  weighs  the  contribution  of  each  pixel  equally.  This 
residual  background  represents  a  small  percentage  of  the  signal  peak,  but  its  effect  on  the 
model  and  measurement  comparison  in  regions  where  the  signal  is  small  can  not  be 
understated.  Subtraction  of  this  constant  background  gives  rise  to  small  regions  where 
the  resulting  intensity  is  negative.  This  unfortunate  effect  is  mitigated  by  setting  all 
negative  values  to  zero. 

Registration  Errors. 

Despite  efforts  to  precisely  align  the  source  image  with  the  center  of  the  focal 
plane,  some  discrepancy  will  inevitably  exist.  To  correct  this  problem,  a  simple  image 
registration  algorithm  is  applied  to  the  model  image.  The  algorithm  shifts  the  model 
image  by  one  pixel  in  each  of  8  separate  directions  and  compares  the  result  to  the 
measured  image  using  the  squared  error  technique.  The  shift  that  provides  the  least 
squared  error  becomes  the  new  center  of  the  model  image  and  the  process  repeats  itself 
until  the  central  pixel  provides  the  least  squared  error. 

There  are  two  potential  drawbacks  to  this  algorithm.  First,  registration  is  limited 
the  precision  provided  by  shifts  made  in  whole  pixel  increments.  Second,  the  algorithm 
has  no  mechanism  for  detennining  whether  the  squared  error  minimum  is  local  or  global 
in  nature.  The  first  issue  can  be  reduced  by  selecting  a  more  sophisticated  algorithm  if 
greater  precision  is  necessary  though  the  results  will  show  that,  at  least  for  this 
incarnation  of  the  model,  this  is  not  the  case.  The  second  issue  is  protected  against  by 
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visual  inspection  of  the  results  (a  significant  misalignment  of  the  images  is  readily 
apparent). 

Prism  Alignment  and  Rotation  Errors. 

Mathematically,  prism  alignment  errors  are  expressed  in  terms  of  the  angle  (p  in 
equation  (3.10)  and  the  surrounding  discussion,  cp  is  treated  as  model  parameter  that  is 
minimized  simultaneously  with  the  rotation  error  angle  using  the  correlation  method 
described  above.  Minimization  of  this  parameter  may  also  suffer  from  the  problem  of 
stumbling  into  local  minima  but,  as  care  is  taken  to  properly  align  the  prisms,  alignment 
errors  of  this  type  are  small. 

Prism  Rotation  angle  errors  are  also  eliminated  parametrically.  As  stated 
previously,  the  circular  symmetry  of  the  source  allows  for  rotation  of  the  composite 
image.  After  an  initial  approximation  of  the  rotation  angle  is  made  in  the  model,  small 
corrections  are  applied  though  an  additional  application  of  the  rotation  routine  discussed 
in  the  previous  chapter.  Like  alignment  error  minimization,  rotation  error  minimization 
is  accomplished  though  the  application  of  the  correlation  technique. 
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VII.  Results,  Analysis,  and  Conclusions 


The  measurement  techniques  and  model  registration  processes  outlined  in  chapter 
VI  are  applied  here  to  four  broadband  (3-5  um)  CT  images.  The  measured  images 
represent  a  400°C  blackbody  with  a  manufacturer  specified  emissivity  of  0.97  or  greater 
(EOI,  2004:  1).  Measurements  were  made  sequentially  over  an  approximate  30  minute 
period  with  identical  camera  settings  (frame  rate,  integration  time,  etcetera)  and 
laboratory  atmospheric  conditions.  All  other  conditions  being  equal,  the  difference 
between  each  image  lies  in  a  change  of  rotation  angle  (90°  each  time)  of  the  prism 
assembly.  Each  measurement  is  presented  along  with  a  corresponding  simulated  image,  a 
model  fit  map,  and  a  brief  explanation  of  the  parameters  selected  to  make  the  fit. 

The  model  fit  map  contains  a  pixel  by  pixel  representation  of  the  agreement 
between  the  simulation  and  the  model.  The  value  reported  for  each  pixel  expresses  the 
probability  of  a  match  in  terms  of  the  number  of  Poisson  standard  deviations  that  the 
model  value  lies  away  from  the  measured  value. 

Results 

Comparison  1:  275°  Rotation  Angle. 

Figures  27,  28,  and  29  contain  (in  order)  the  CT  measurement,  the  simulation,  and 
the  fit  map.  The  best  fit  was  arrived  at  using  a  background  bias  of  320  counts  (2.2%  of 
the  peak)  and  a  prism  alignment  error  angle  of  0.1°  (the  rotation  angle  was  also  fit 
parametrically).  Registration  required  a  net  shift  of  2  pixels.  Overall  correlation  between 
the  two  images  is  0.9891. 
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Figure  27.  275°  Prism  Rotation  Angle  Data 
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Figure  28.  275°  Prism  Rotation  Angle  Simulation 
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Figure  29.  275°  Prism  Rotation  Angle  Fit  Map 

The  peak  error  shown  in  figure  29  is  8.53  standard  deviations  and  the  average 
error  in  the  region  specifically  containing  the  signal  is  1.87  standard  deviations. 

Comparison  2: 184.5°  Rotation  Angle. 

Figures  30,  31,  and  32  contain  the  CT  measurement,  the  simulation,  and  the  fit 
map  of  the  second  set  of  measurements.  The  best  fit  was  arrived  at  using  a  background 
bias  of  390  counts  (3.02%  of  the  peak)  and  a  prism  alignment  error  angle  of  0°. 
Registration  required  a  net  shift  of  6  pixels.  Overall  correlation  between  the  two  images 
is  0.9894. 
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Figure  30.  184.5°  Prism  Rotation  Angle  Data 
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Figure  31.  1 84.5°  Prism  Rotation  Angle  Model 
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Figure  32.  184.5°  Prism  Rotation  Angle  Fit  Map 
The  peak  error  shown  in  figure  32  is  7.15  standard  deviations  and  the  average  error  in  the 
region  specifically  containing  the  signal  is  1.35  standard  deviations. 

Comparison  3: 100°  Rotation  Angle. 

Figures  33,  34,  and  35  contain  the  CT  measurement,  the  simulation,  and  the  fit 
map  of  the  third  set  of  measurements.  The  best  fit  was  arrived  at  using  a  background  bias 
of  380  counts  (3.00%  of  the  peak)  and  a  prism  alignment  error  angle  of  0.1°. 

Registration  required  a  net  shift  of  13  pixels.  Overall  correlation  between  the  two  images 
is  0.9895. 
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Figure  34.  100°  Prism  Rotation  Angle  Simulation 


Figure  33.  100  Prism  Rotation  Angle  Data 
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Figure  35.  100°  Prism  Rotation  Angle  Fit  Map 

The  peak  error  shown  in  figure  35  is  7.1 1  standard  deviations  and  the  average  error  in  the 
region  specifically  containing  the  signal  is  1 .65  standard  deviations. 

Comparison  4: 14.5°  Rotation  Angle. 

Figures  36,  37,  and  38  contain  the  CT  measurement,  the  simulation,  and  the  fit 
map  of  the  fourth  set  of  measurements.  The  best  fit  was  arrived  at  using  a  background 
bias  of  340  counts  (2.44%  of  the  peak)  and  a  prism  alignment  error  angle  of  0.2°. 
Registration  required  a  net  shift  of  3  pixels.  Overall  correlation  between  the  two  images 
is  0.9901. 
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Figure  37.  14.5°  Prism  Rotation  Angle  Simulation 


Figure  36.  14.5  Prism  Rotation  Angle  Data 
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Figure  38.  14.5°  Prism  Rotation  Angle  Fit  Map 

The  peak  error  shown  in  figure  38  is  7.62  standard  deviations  and  the  average  error  in  the 
region  specifically  containing  the  signal  is  1.76  standard  deviations. 

Summary  of  the  Experimental  Data. 

Table  3  contains  a  consolidated  version  of  the  information  supplied  above  to 
facilitate  the  discussion  in  the  following  sections. 

Table  3.  Results  from  Four  Comparisons 


Comparison 

Bias 

(counts) 

Rotation 

Angle 

(deg) 

Alignment 

Error 

(deg) 

Correlation 

Coefficient 

Max 

Error 

(std 

dev) 

Mean 

Error 

(std 

dev) 

1 

320 

275.0 

0.1 

0.9891 

8.53 

1.87 

2 

390 

184.5 

0.0 

0.9894 

7.15 

1.35 

3 

380 

100.0 

0.1 

0.9895 

7.11 

1.65 

4 

340 

14.5 

0.2 

0.9901 

7.62 

1.76 
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Analysis 

With  the  exception  of  rotation  angle,  the  comparison  results  are  similar  for  each 
of  the  four  data  sets.  The  following  section  describes  how  each  of  the  columns  in  table  3 
can  be  used  to  reinforce  or  detract  from  this  opinion.  Conclusions  based  on  this  analysis 
are  presented  in  the  following  section. 

The  background  bias  floats  between  320  and  390  counts  or  2.20%  to  3.02%  of  the 
signal  peak.  The  relative  size  and  consistency  in  background  bias  helps  to  reaffirm  the 
initial  justification  for  including  this  parameter,  which  was  to  minimize  the  impact  of 
emission  from  the  area  immediately  surrounding  the  blackbody  cavity.  Recalling  that  the 
standard  deviation  of  the  Poisson  distribution  goes  with  square  root  of  the  measured 
value,  an  extra  300  counts  could  have  easily  been  interpreted  as  additional  error  in  places 
where  no  such  error  existed. 

Surprisingly,  the  only  apparent  major  discrepancy  in  the  experiment  comes  from 
rotation  angle.  Following  the  order  in  which  the  measurements  are  presented,  each 
subsequent  measurement  should  represent  a  net  rotation  of  90°.  Rotation  between  the 
first  two  measurements  is  almost  exactly  that  but  the  angular  difference  between  the 
second  and  third  measurements  is  85.5°.  Confounding  matters  further,  the  difference 
between  the  third  and  fourth  measurements  is  1 14.5°.  This  discrepancy  is  apparently  a 
motor  control  issue  and,  since  the  model  can  compensate  for  any  rotation  angle,  its 
impact  on  this  research  is  negligible. 

Prism  alignment  error,  the  apparent  angular  difference  between  the  theoretical 
orientation  of  each  prism  and  the  actual  orientation,  varies  slightly  across  the  four 
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measurements  (0  to  0.2  degrees).  Since  no  adjustments  were  made  to  the  prism  assembly 
in  the  interval  between  measurements,  two  possible  explanations  for  this  shift  exist. 

First,  a  physical  shift  may  occur  during  prism  rotation  though  this  seems  unlikely  as  both 
prisms  are  held  in  place  via  set  screws.  Second,  the  shift  represents  numerical  error  in  the 
model  and  measurement  reconciliation  software.  Given  that  the  relative  error  is  small 
and  that  it  appears  to  follow  no  particular  pattern,  the  second  explanation  would  appear  to 
be  the  most  plausible. 

Though  all  four  correlations  are  quite  high,  the  correlation  between  the  measured 
and  simulated  images  do  not  speak  to  the  overall  merit  of  the  simulation  based  upon 
magnitude  alone  (anecdotal  evidence  of  this  statement  is  provided  in  the  example  given 
in  the  previous  chapter).  Rather,  the  correlation’s  relevance  to  this  research  is  derived 
from  its  consistency  across  the  four  measurements.  Note  from  table  3  that  the  maximum 
difference  between  any  two  correlations  is  limited  to  0.001.  By  itself,  this  piece  of 
information  is  an  indication  of  consistency  in  the  results,  but  not  necessarily  an  indication 
of  accuracy. 

Two  important  points  come  out  of  the  measurement  of  maximum  error  across  the 
four  sets  of  results.  First,  the  max  error  is  spread  out  between  7.11  and  8.53  standard 
deviations,  which  is  another  indication  of  consistency  throughout  the  results.  Second,  the 
max  error  is  always  found  in  the  same  regions  of  dispersion  as  seen  in  figures  29,  31,  34, 
and  37.  In  other  words,  the  maximum  error  regions  (represented  by  two  bright  lobes  in 
the  figures)  correspond  to  regions  of  the  focal  plane  affected  by  dispersion  in  the  same 
spectral  region. 


73 


This  error  may  be  the  result  of  a  variation  in  focal  length  of  the  three  lens  system 
as  a  function  of  wavelength.  This  variation  will  manifest  itself  physically  as  a  gradient  in 
magnification.  A  return  to  geometric  optics  helps  to  test  this  hypothesis.  Focal  length,  f, 


is  detennined  for  thin  lenses  in  air  using  the  lens  maker’s  equation  (Hecht,  2002: 158) 


7  =  fa"1) 


VRi  R2y 


(7.1) 


where  Rx  represents  the  curvature  of  the  lens  and  ni  is  the  refractive  index  of  the  lens. 

The  in  focus  magnification  from  equation  (3.28)  is  approximated  as  a  function  of 
wavelength  under  the  assumption  that  the  manufacturer’s  specified  focal  length  occurs  at 
4  um  (the  center  of  the  CT  bandpass).  Figure  39  is  a  plot  of  approximate  magnification 
as  a  function  of  wavelength. 
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Figure  39.  Magnification  as  a  Function  of  Wavelength 

The  overall  change  in  magnification  is  subtle  but  certainly  enough  to  affect  the 
outcome  of  the  simulation.  Note  that  figure  39  approximates  the  in  focus  magnification 
over  the  bandpass;  in  reality,  the  sensor  can  only  be  truly  in  focus  at  one  wavelength. 
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Consequently,  figure  39  is  useful  for  demonstrating  that  the  effect  exists,  but  provides 
nothing  practical  to  add  to  the  simulation. 


Alternatively,  the  apparent  change  in  magnification  may  be  the  result  of 
distortion,  which  is  defined  to  be  a  variation  in  transverse  magnification  with  respect  to 
off  axis  image  position  (Hecht,  2002:  266).  Though  it  can’t  be  ruled  out,  the  distortion 
argument  is  more  difficult  to  justify  because  the  regions  of  greatest  error  appear  to 
correspond  with  a  change  in  spectral  rather  than  spatial  positioning. 

Mean  error  is  also  consistent  throughout  the  series  of  comparisons.  This  figure  of 
merit  is  bracketed  between  1.35  and  1.87  standard  deviations,  though  this  range  may  be 
somewhat  artificially  depressed.  Inspection  of  each  fit  map  shows  apparent  curves  of 
near  zero  error  that  wind  throughout  the  scene.  Though  the  pixels  values  of  the 
simulation  and  measurements  are  nearly  equal  in  these  regions,  this  is  probably  due  to 
coincidental  overlaps  in  the  two  intensity  patterns  rather  than  the  manifestations  of  an 
exact  model.  Figure  38  helps  to  illustrate  this  concept. 
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The  two  imaginary  intensity  profiles  in  figure  38  overlap  in  two  locations.  This  overlap 
creates  regions  where  the  residual  is  numerically  small  even  though  the  overall  fit  is  no 
better  in  these  regions  than  in  anywhere  else.  This  effect,  extrapolated  into  three 
dimensions,  is  readily  apparent  in  all  four  of  the  fit  maps  and  is  responsible  for  lowering 
the  average  error  though  its  overall  impact  is  difficult  to  assess. 

Conclusions 

Research  Summary. 

The  purpose  of  this  research  was  to  demonstrate  the  validity  of  a  phase 
screen  modeling  approach  by  simulating  a  chromotomographic  hyperspectral  imaging 
system.  The  progression  of  this  document  mirrors  the  chronological  progression  of  the 
research:  constructing  the  CT  sensor  (chapters  II  and  III),  modeling  the  sensor  (chapters 
IV  and  V),  and  making  measurements  for  comparison  (chapters  VI  and  VII).  Ultimately, 
success  or  failure  of  the  model  was  based  on  a  measure  of  its  statistical  agreement  with 
the  measured  data. 

Assumptions. 

Three  primary  assumptions  were  built  into  the  modeling  process.  The  first  (and 
predominant)  assumption  was  that  the  CT  system  is  linear  and  invariant  in  intensity  over 
the  region  of  the  focal  plane  occupied  by  the  image.  This  assumption  is  critical  if  the 
Fourier  transform  model  of  propagation  is  valid.  Second,  radiation  entering  the  focusing 
lens  can  be  modeled  as  a  corruption  of  a  collimated  beam.  This  assumption,  which 
primarily  ignores  the  effects  of  diffraction  during  the  collimation  process,  paved  the  way 
for  the  single  lens  approximation.  Third,  all  aberrations,  whether  they  be  from  alignment 
of  or  innate  to  the  optical  system  or  its  components,  can  be  captured  and  reproduced  in 
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the  phase  screen  measurements.  In  one  sense,  this  assumption  is  a  restatement  of  the  first 
but  it  extends  the  definition  to  include  radiation  of  all  wavelengths. 

In  addition  to  the  primary  assumptions,  which  are  unavoidable,  several  secondary 
assumptions  crept  into  the  research.  These  assumptions  include: 

(1)  Any  deviation  from  perfect  prism  alignment  is  described  in  the  model. 

(2)  The  blackbody  source  is  uniform  and  Lambertian. 

(3)  All  photons  entering  the  leading  aperture  of  the  sensor  are  accounted  for. 

(4)  Reflected  radiation  can  be  ignored 

Conclusions. 

This  research  has  shown  that,  over  some  field  of  view  and  some  spectral  range, 
the  phase  screens  generated  using  the  phase  retrieval  process  can  be  used  to  approximate 
the  intensity  patterns  produced  by  the  chromotomographic  imaging  system.  The  overall 
agreement  between  the  measured  and  simulated  data  (expressed  in  tenns  of  the  generally 
low  mean  error)  is  evidence  of  this  statement. 

Regions  where  the  simulated  and  measured  data  disagree  most  appear  to  be 
dominated  by  either  an  off-axis  magnification  effect  (distortion),  or  a  chromatic  variation 
in  magnification.  These  potential  sources  of  error  can  be  reduced  or  eliminated  by 
adjusting  the  sensor  (achromatic  lenses,  shifting  or  reducing  the  field  stop)  or  by 
generating  multiple  phase  screens  at  regular  intervals  to  compensate  for  these  effects. 

The  successful  assimilation  of  the  phase  screen  technique  into  the  Fourier 
propagation  model  provides  a  conceptual  basis  for  applying  the  same  technique  to  a 
modified  CT  image  reconstruction  algorithm. 
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Recommendations  for  Future  Research. 

The  conclusions  statement,  statement  of  assumptions,  and  original  motivation  for 
this  research  all  lend  themselves  immediately  to  a  list  of  future  research  projects.  If  the 
phase  screens  are  to  be  imported  into  a  chromotomographic  deconvolution  algorithm, 
then  the  next  appropriate  research  step  is  to  identify  the  largest  field  of  view  over  which 
the  primary  assumptions  apply.  Concurrent  efforts  should  be  spent  identifying  the  useful 
spectral  bandwidth  of  any  particular  phase  screen. 

Recall  that  the  Gerchberg-Saxton  phase  retrieval  algorithm  could  not  generate  a 
phase  screen  that  perfectly  reproduces  the  Richardson-Lucy  point  spread  function.  While 
the  utility  of  Gerchberg-Saxton  is  well  known,  other  algorithms  exist  that  can  perform  the 
same  function  (Schulz,  1992:  1266).  Future  research  should  concentrate  on  selecting 
another  phase  retrieval  algorithm  or  modifying  an  existing  one  to  suit  the  specific  needs 
of  the  CT  optical  system. 

Modeling  research  could  continue  by  expanding  the  single  propagation  model  into 
a  multi-propagation  model.  This  step  would  increase  the  complexity  of  the  phase 
retrieval  process  immensely  by  requiring  that  a  separate  phase  screen  be  provided  for 
each  element  in  the  optical  system.  Since  the  phase  screen  conveys  both  innate  and 
alignment  aberrations,  each  screen  after  the  focusing  lens  would  have  to  be  detennined 
by  some  sort  of  deductive  process.  If  successful,  this  research  step  would  effectively 
eliminate  the  second  primary  assumption. 
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